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ABSTRACT 


An important branch of mathematical programming 1s concerned with 
optimization in systems described by networks. This paper describes an 
integrated suite of advanced techniques for dealing with minimum cost 
network flow formulations. Written in Pascal and implemented on a 
microcomputer representative of current small computer technology (the 
APPLE II), this package places unprecedented modeling versatility and 
solution capadility on the analyst's desktop. Able to solve small to 
medium size problems (3000 arcs or less) at reasonable speeds, programs 
to handle capacitated linear, nonlinear (convex separable), mixed 
integer and elastic ranged linear models in addition to comprehensive 
contro! and data management routines are included. Problem size and 


solution speed benchmarks are given for a variety of models. 
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NOTATION 


Except as otherwise indicated in the text, the following notational 


conventions have been used for all mathematical expressions: 


Scalars Lower 


case Greek and Latin letters (e.g., a, a). 


Vectors Lower case Latin letters with a bar above (a). 

Vector Components Lower case Latin letters with a lower case 
Latin subscript (a.). 

Matrices Upper case Latin letters (A). 

Matrix element Lower case Latin letters with lower case 
Latin subscripts (a, 5). 

Transpose Superscript T (aly. 

Sequences Enclosed in braces ({ay}). 
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I. INTRODUCTION 


Mathematical programming can be defined as the use of mathematical 
representations (models) to plan (program) an allocation of scarce 
resources among competing activities [Ref. 1]. An important branch 
of mathematical ae ne deals with optimization in systems described 
by networks or collections of points (nodes) connected by links (arcs). 
Such models arise explicitly in a variety of applications and include 
many familiar distribution and transportation problems. Bradley [Ref. 2] 
suggests that network models have been so widely used because: (1) they 
accurately model many applications, (2) they are more readily accepted by 
nonanalysts than other models (they pictorially resemble the physical 
process being modeled), and (3) efficient algorithms are available. 
Additionally, many models can easily be transformed into equivalent 
network representations by direct manipulation (e.g., assignment problems) 
or exploitation of a primal-dual relationship (e.g., critical path 
problems). The object of all these formulations is usually to minimize 
the cost of moving a single commodity through the network. The general 
form of this problem, the network programming problem (NPP), is 

(NPP)  MINIMIZE f(x) (cost function) 


Se te Ax = b (node flow-conservation constraints) 


1<x<u_ (are flow bounds), 
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where x ig an n-dimensional vector 
(n is the number of arcs), 
x. represents the magnitude of the flow on arc i, 
A is anm by n (m is the number of nodes) 
matrix defined as follows: 
+] if arc j is directed away from node i 
aij ~ -l1 if arc j is directed toward node 1 
QO otherwise, 
b is an m-dimensional vector of flow requirements at 


each node, and 


f(x) is some function that relates cost to arc flow. 


Although NPP can be solved by classical general purpose constrained 
methods (depending on the exact form of f(x)), better techniques exist. 
Specialized data structures employed in conjuction with modifications 
of traditional optimization procedures have resulted in the development 
of extremely efficient algorithms for the solution of NPP [e.g., Ref. 3]. 
This work has been motivated in part by the intrinsic usefulness and wide 
applicability of these models and further stimulated by the increased 
availability of computational devices. Due to these advances, the 
operations analyst is now able to represent larger network models, and 
answer questions concerning their optimal flow easier than ever before. 

As interest in network programming continues to grow, the supply 
Of applications software will undoubtedly keep pace. This has been the 
case to date: however, the vast majority of the emerging software appears 
very specific in nature and tailored to particular classes of network 
models. Certainly there are many superb codes available, but if, for 
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example, the capability to work with both linear and nonlinear cost 
functions is desired, then the employment of two separate and possibly 
quite different programs is required. Even a relatively mundane task 
such as linearizing a nonlinear model to obtain a feasible starting 
solution point either requires incorporation of extra code or an exercise 
in data manipulation and program linking. Such shortcomings indicate 
that unified network flow programming packages able to cope with several 
different types of models would be quite useful. 

Concurrent with these advances in network programming, computer 
technology has experienced breath-taking progress. Not only have com- 
puters improved in sheer power, but they have become increasingly 
compact and less expensive (relative to capability). A dramatic example 
of this data-processing revolution is the evolution of desk-top computers. 
These small computers, the so-called microcomputers, appeared on the 
scene in the early 1970's and today can provide computer power comparable 
to some room-sized machines of a decade ago at a tiny fraction of the 
cost. More importantly, these devices are becoming so common that their 
accessibility to technical personnel is forecasted to be virtually 
universal in the near future [Ref. 4]. E£. M. L. Beale [Ref. 5] commented 
on this dissemination of computer power and predicted increased use of 
microcomputers by operations analysts to solve smaller models locally 
with reliance on computer centers for large projects. Tanenbaum [Ref. 6] 
suggests that the cost of small computers relative to communication 
expenses now makes it attractive to analyze data at its source and send 


only summaries back to large computers via networking arrangements. 
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Even though the hardware is available and its potential clearly 
recognized, suprisingly little operations analysis software of reasonable 
complexity and generality has been reported for microcomputers. Aside 
from a few decision analysis programs [e.g., Ref. 7, 8, 9], effort has 
been mainly concentrated on statistical applications. Morgenson 
[Ref. 10] and Isbell [Ref. ll], for example, have developed extensive 
data analysis packages incorporating sophisticated techniques. Addition- 
ally, there are numerous commercial statistical products of varying 
quality and capability for microcomputers. 

Surprisingly, microcomputer-based mathematical programming software 
is still virtually non-existant, despite the obvious utility of such 
programs. Undoubtedly, there are many ad hoc, rudimentary implementa- 
tions of basic methods; however, no reference to any work of consequence 
has been found in the open literature. Feasibility studies [Ref. 12; 13] 
verify the desirability of optimization packages for small computers, but 
only explore rather elementary network algorithms. These studies further 
suggest that network flow problems are likely candidates for microcomputer 
solution because of the efficient algorithms at hand. Economic justifi- 
cation for employing the microcomputer for optimization purposes has not 
been established conclusively. One attempt to favorably compare such use 
with the alternative of large, general purpose computer systems [Ref. 13] 
has been severely criticized [Ref. 14] and the issue remains undecided. 

The research reported here investigates the construction of unified 
network flow optimization packages, the mathematical programming potential 


of microcemputers, and (secondarily) the economic feasibility of such 


v/ 





devices with respect to optimization applications. This work was under- 
taken with the following (superficially disparate) goals: 

1. Development of a unified, versatile network-flow optimization 
package capable of handling a variety of models and utilizing 
"State-of-the-art" algorithms. 

2. Implementation of this package on a widely available micro- 
computer to explore the usefulness of smaller computers in an 
Optimization context. 

Portions of this work were presented at the CORS/ORSA/TIMS joint 

meeting in Toronto, May, 1981. At that international meeting, attendees 
expressed surprise that such an integrated network optimization package 


used a microcomputer as a host. 
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IT. DESIGN AND IMPLEMENTATION CONSIDERATIONS 


Microcomputer implementation of large programming projects requires 
careful consideration of every aspect of program design in order to fully 
exploit the limited resources available. This section briefly discusses 
the design criteria utilized and outlines the rationale for the decisions 


that shaped the package. 


A. BACKGROUND 

Invariably many different program configurations can be constructed 
to accomplish a given task, however, some approaches are better than 
others. Program effectiveness is a well-studied field and the literature 
contains many characterizations of superior software designs. Kreitzberg 
and Shneiderman [Ref. 15], for instance, define a "good" program as one 
that is correct (provides desired results), fast, accurate, hardware inde- 
pendent, efficient with storage, and easily modified. Thenson [Ref. 16] 
and Lientz [Ref. 17] further suggest that operations analysis software 
should place heavy emphasis on "user impact" or ease of use. These 
attributes, commendable as they may be, are of little value if the 
program cannot solve meaningful problems in a reasonable amount of time. 
Meaningful, in this context, includes not only the intrinsic usefulness 
and generality of the modeling facilities provided, but also the size and 
complexity of the representable formulations. Since it is highly unlikely 


that a procram can be optimal in every respect, compromises are inevitable. 


IS, 








Numerous factors interact to establish program limitations with 
the inherent capabilities of the host computer playing the dominant 
role. Choice of algorithms, data structures, programming language, 
memory and peripheral management, and programming tactics are also 
Significant determinants of ultimate software efficiency. For the 
traditional user of large computer systems, the improper choice of one or 
more of these variables usually results in merely a slightly degraded 
package rather than outright failure. The small-computer environment is 
less forgiving and does not allow this luxury without severe performance 


penalties. 


B. HARDWARE 
1. Microcomputer Characteristics 

A typical microcomputer system consists of a central processing 
unit, memory, peripheral devices, and software. Exclusive of peripherals, 
such packages physically resemble an electric typewriter. The specific 
form of the processing unit varies, but most can directly address 65,536 
eight-bit words of random-access memory. This “fast" (typically 200-400 
nsec) or core memory is usually supplemented by slower (but more plentiful) 
storage (worst case, 1-2 seconds seek time) in the form of disk drives. Each 
disk drive provides approximately one hundred thousand (or more) words of 
memory space on removable magnetic media. The number of disk drives allowed, 
their individual capacities, and access times depend on the microcomputer. 
More sophisticated (and expensive) offline memory devices are also available, 
capable of storing millions of words. Common peripherals include communica- 


tion equipment to access other computers, video displays, and printers. 
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2. layget Microcomputer 


Clearly, as technology is changing so rapidly, any attempt to 
identify an optimal hardware selection would be futile. Instead, an 
acceptable system representative of available products, the APPLE it 
microcomputer, 1S considered. All software development has been per- 


formed on an APPLE II system configured as shown in Figure l. 


1. APPLE II microcomputer with 65,000 words 
of memory. 


2. Two disk-drive offline storage devices 
(150,000 words each). 

3. UCSD Pascal (language card). 

4. Communications device (modem). 

5. Printer. 

6. Eighty character by twenty-four line 


video display. 


Fig. 1. Development Hardware 


Although this machine is not the most powerful of its class, it 
is a reasonable choice for many operation analysis applications due to 
its availability, capability, and low cost. Introduced in 1977, the 
APPLE II is a widely-used device, still in production, and likely to 
remain so for the next few years [Ref. 18]. As over 200,000 units have 
been manufactured and distributed worldwide, a broad base of technical 
Support exists and numerous firms offer peripheral products. The APPLE 
II can support all of the popular microcomputer programming languages 
(BASIC, Pascal, FORTRAN, and Assembly language). Easily expandable, 


LAPPLE II is a trademark of Apple Computer Inc. 
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complete systems can be purchased for between two and four thousand 
dollars. 
3. Machine Independence 

All computers, regardless of size, have unique features and one 
must be careful if machine independence of software is desired. Given the 
variety of microcomputers available today and the rapid changes expected, 
portability between existing and conceptual machines is necessary if the 
considerable investment involved in the development of sophisticated 
programs is to be protected. Thus the importance of machine independence 
should not be underestimated. For this reason, features peculiar to the 
APPLE II have been avoided wherever possible. Thus the software presented 
should run, with minimal modifications, on any equivalent hardware 


package. 


C. MODEL AND SOLUTION TECHNIQUE SELECTION 

The determination of specific processes to model and subsequent 
selection of solution techniques compatible with computing resources are 
critical issues in the design of any optimization software. 

1. Models 

Once the general class of problems has been established, a few 

carefully chosen models often can adequately represent the majority of 
anticipated situations. Specialized models can then be added as necessary. 
The following minimum cost network flow models are considered essential 
for any comprehensive package: 

Oo Linear cost function with bounded variables. 


o Nonlinear convex separable cost function with bounded variables. 
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As useful as these two models are, they cannot depict some 
desired formulations. Additional flexibility is provided by: 
Oo Linear cost function with bounded variables and elastic ranged 
constraints (see Section III, Subsection C). 
o Linear cost function with elastic range constraints and bounded 
variables, any of which may also be specified as “l-u" variables 


(only allowable flow is at one of the bounds). 


These four models comprise a basic package capable of repre- 
senting a wide variety of single commodity network flow problems. 
Included are fundamental examples of linear, nonlinear, and mixed-integer 
network optimization models. Building from this network paradigm, these 
features can be used to represent myriad mathematical models. Other 
models could easily be added, however, time constraints precluded further 
extension of the package. 

2. Solution Methods 

When dealing with microcomputers, algorithms are difficult 
to select. Simple algorithms require little storage, but are typically 
inefficient and slow. More complicated approaches offer speed at the 
expense of increased memory usage and expanded data structure requirements. 
Large computers are able to use the efficient algorithms because the 
additional storage requirements are infinitesimal compared to the total 
memory available. For small machines, the choice is not so clear. A 
large portion of memory is consumed by the added complexity of the 
advanced algorithms, resulting in significant reduction of maximum problem 
size. If one chooses the simpler methods, larger programs are accommo- 
dated but processing time can soar to unacceptable levels for all problems. 
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Execution speed appears to be far more important than maximum 
problem size. Providing quick answers to small problems is the most 
likely optimization role for microcomputers at the present time. Very 
large problems will continue to require extensive computational facilities 
well beyond the capabilities of current microcomputers. Given the 
predicted advances in small computers, memory-saving techniques should 
decrease in importance as devices with greater storage capacities become 
available. Using inefficient methods instead of the more capable but 
memory-intensive algorithms does not seem to be justified. 

Numerical representation must also be considered. Small computers 
are usually more restricted than their larger counterparts in (1) mantissa 
precision, (2) real variable exponent range, and (3) maximum allowable 
integer size. If provisions exist (or can be created) to deal with these 
problems, performance degradation may result. Thus candidate algorithms 
for microcomputer use must be relatively insensitive to such limitations. 

The programs selected for these models are considered state-of- 
the-art and all represent current research efforts. Utilizing advanced 
algorithms and efficient data structures, these programs have been 
included on the basis of their suitability for microcomputer adaptation 
and mutual compatability. The original programs were coded in FORTRAN 
for implementation on large computers. Each program has been translated 
into UCSD Pascal and modified as necessary to enhance its efficiency in a 
microcomputer environment. 

GNET (Ref. 3] is the foundation of the package. A highly regarded 
and widely-used code, it is one of the fastest methods available for 


solving linear minimum cost network models with bounded variables. 
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Additionally, GNET can be used to linearize nonlinear problems and 
produce initial feasible solutions. The documented speed of the method 
[Ref. 2] and its efficient data structures make it ideal for microcomputer 
implementation. 

Nonlinear formulations are solved by NLPNET [Ref. 19], a program 
by Dembo that utilizes data structures in the spirit of GNET and new 
adaptable direction-finding techniques. The user of NLPNET can control 
the accuracy of the direction-finding process and thus the amount of 
computation effort expended for a solution. This feature makes such an 
algorithm perfect for microcomputers. 

Elastic ranged constraints are handled by ENET, a recently 
introduced code by Brown and Graves [Ref. 20]. An extension of GNET, it 
uses the same efficient methods and data structures. 

Finally, UNET, another new code by Brown and Graves [Ref. 20], is 
used to solve elastic ranged constraint problems with bounded variables, 
some of which may be specified as "l-u" variables (flow at either bound, 
only). Also employing data structures similar to GNET, this algorithm 
requires successive calls to ENET to solve enumeration subproblems. 

Although ENET can duplicate the capabilities of GNET, the latter 
is explicitly retained to solve larger problems of the simpler structure 
amenable to GNET. This is not an algorithmic disadvantage of ENET, but a 
consequence of the package design. Since UNET repeatedly calls ENET, 
both procedures are kept in memory together to reduce disk access time. 
Therefore, as the programs are constructed, GNET 1s a more compact code 
than the ENET/UNET combination. It is intended that future revisions of 


the package will employ only ENET and UNET. 
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D. PROGRAMMING LANGUAGE 

Pascal, specifically the University of California at San Diego (UCSD) 
version, appears to currently offer the best programming language facil- 
ities for large-scale microcomputer programming projects. Possible 
alternatives are BASIC and FORTRAN, languages that suffer from crippling 
Shortcomings as presently implemented on microcomputers. 

1. BASIC 

Once the only high-level language available for small computers, 

BASIC remains very popular. While BASIC is perhaps suitable for some 
applications, serious mathematical programming efforts are hampered by 
its limited speed, primitive programming power, and lack of transport- 
ability. When installed as an interpreted language (the usual case), 
each source statement must be translated into machine language every time 
it is encountered in the logical flow of the program. This results in 
long execution times, especially for complex programs with many iterative 
Structures (i.e., optimization programs). Additionally, most dialects of 
BASIC allow only global (accessable to all program portions) variables 
with restrictive naming conventions. Variable name conflicts often arise 
in such circumstances greatly reducing the portability of subprograms. 
The use of subprograms is further hindered by the inability to pass 
parameters between program fragments without resorting to global variable 
reassignments. Compiled versions (one-time translations of source 
Statements to machine language) of BASIC are available that supposedly 
alleviate these problems, however, they are necessarily machine-specific, 
further aggravating the transportability issue. There are probably as 


many different versions of BASIC on the market as there are types of 
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microcomputers. Hardware manufacturers modify the language to fit their 
particular needs resulting in potential difficulties when programs are 
moved to different machines. It 1s indeed a rare BASIC program that can 
be transported without modification, with the chances of such success 
diminishing as program length and complexity increases. 
2. FORTRAN 
Implementations of FORTRAN for microcomputers, although increasing 
in number, are still rather uncommon. A compiled language, FORTRAN 
compares favorably to other languages in terms of speed and programming 
power. The biggest drawback of FORTRAN is lack of standardization, a 
hinderance that has plagued FORTRAN users on large computers for years. 
Compiled code is completely machine-specific and thus not portable. 
Transportability of the source statements varies depending of the degree 
of similarity between the two versions of FORTRAN in question. This lack 
of standardization combined with the inability to exchange compiled code 
between different microcomputers makes FORTRAN a poor choice if program 
portability is a goal. 
3. Pascal 
a. Standard Pascal 

Pascal is a relatively new language (compared to BASIC 
and FORTRAN) having been formally defined by N. Wirth in 1971 [Ref. 21]. 
Named after the famous mathematician Blaise Pascal, the language was 
Originally intended as a vehicle to teach computer programming. An 
excellent instructional tool, Pascal is now recognized as a powerful 
general-purpose language. Similar to ALGOL, Pascal provides a rich 


set of program and data structuring tools. Thus a great portion of 
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the housekeeping tedium of programming 1s assumed by the language itself 
(in the form of the compilation process) and is totally transparent to 
the user. Some of these facilities can be simulated in other languages 
through skillful programming, but this adds complexity to the development 
process. Pascal, by admitting long variable names and free-form coding, 
is a self-documenting language further reducing development and mainte- 
nance effort. Totally modular in concept, programs consist of blocks of 
similarly structured code with true local and global variables allowed at 
all levels. A Pascal program is therefore a logical collection of 
independent modules that allows direct translation of complicated algo- 
rithms into source code by partitioning the problem into smaller and more 
manageable pieces. This divide and conquer technique, commonly known as 
"structured programming," is an effective method to address large projects. 
Figure 2 is an example of this block structure. Extension of Pascal 
programs is a Simple matter of adding modules of code and linking them to 
the original program via parameters. Pascal also requires explicit type 
definition of all variables, thus minimizing the chance of a non-local 


typographical mistake introducing a program error. 


Block Heading (Program, procedure, or function) 
Variable declarations 
Local blocks 


Begin body of block 
Executable statements 
End body of block 


Fig. 2. Pascal Block Structure 
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b. UCSD Pascal 
These impressive programming features alone are enough to 
warrant the use of Pascal for difficult microcomputer programming projects. 
When one examines the UCSD version, however, the evidence seems over- 


whelming. An extension of standard Pascal, UCSD Pascal® 


was developed 

by a group headed by K. L. Bowles in 1975. Explicitly designed to run in 
a small computer environment, it retains most of the features of standard 
Pascal while adding many extensions that increase the power of the 
language. For our applications, we are interested in two of these 
additional facilities that (1) enhance memory manipulation capabilities 
and (2) provide for extremely transportable programs. 

The memory management techniques offered by UCSD Pascal are 
truly remarkable for a microcomputer software system. The programmer has 
complete control over which portions of code are to be in memory at a 
given point in the program execution. For small programs, the entire 
code is generally loaded into memory at the beginning of execution while 
larger programs can be segmented to provide the desired partitions. A 
small portion of the program, again totally definable by the programmer, 
must be in residence at all times to control the overlaying process. 
Complete programs can be called into memory, each with its own local 
memory control, executed, and other programs subsequently activated. In 
this manner, programs can be chained together to form coherent packages. 


Additionally, groups of often-used blocks of code can be pre-compiled and 


2UCSD Pascal is a trademark of the Regents of the University of 
California. 





29 





placed in libraries. When activated, the specific library code is loaded 
into memory with the calling program and executed. Library routines can 
either remain in memory for the entire program execution process or be 
discarded after each use (at the programmer's option). These methods 
allow programs that would otherwise be impossible to implement on a 
microcomputer--the package presented is such a case as the combined 

code of the component programs far exceeds the memory capacity of most 
microcomputers. 

UCSD Pascal is a hybrid compiled/interpreted language, a 
concept that produces unprecedented machine independence. Source state- 
ments are first compiled to an intermediate pseudo (p) code that is then 
executed by a machine-dependent interpreter. This greatly enhances 
portability as the interpreter, which normally takes about six man-months 
to write [Ref. 22], is the only part of the system that must be changed 
to take advantage of a new hardware configuration. The original compiled 
source code can be transported without modification to any machine that 
has a p-code interpreter installed. Figure 3 depicts the relationship 
between the components of the p-code system (Lewis [Ref. 23]). Such 
interpreters now exist for all the commonly used microcomputer central 
processing units [Ref. 24]. The use of p-code is not restricted to UCSD 
Pascal; other languages (e.g., FORTRAN-77 and PL/I) have been treated in 
this manner to achieve similar portability. However, UCSD Pascal is 
currently the only p-code implementation that is widely available. 

As previously indicated, speed is an important programming 


language attribute. The hybrid nature of UCSD Pascal results in 
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Fig. 3. Relationship Between P-Code and Pascal 


execution time somewhat slower than pure-compiled languages such as FORTRAN, 
but much faster than interpreted languages. Gagne [Ref. 25] reports that 


Pascal runs three times as fast as the very best BASIC and ten to thirty 





times faster than most. Informal timing tests comparing Pascal with APPLE II 

BASIC confirms these performance assertions. Any loss of speed through the 
interpretation of p-code is more than offset by the other desirable features 
of the language, especially portability. 

c. Limitations of Pascal 
Pascal is not without its difficulties--some are minor 

irritants while others are more important. Foremost is the precision of 
real variables allowed by the APPLE II version of UCSD Pascal. Without 
taking any special programming measures, six to seven significant figures 


can be accommodated in a mantissa representation equivalent to IBM 360/370 


Silk 








Single precision. It must be understood that this is an implementation 
issue on the APPLE II and not an intrinsic limitation of Pascal. By 
suitable programming, the precision can be extended; however, it was 
decided not to do so because of the speed degradation that would surely 
have resulted. The only package module that utilizes vast amounts of 
real arithmetic is NLPNET, the nonlinear network code. As shall be 
shown, for the nonlinear problems that have been attempted, limited 
precision does not appear to hinder performance. 

Another disturbing omission is that all arrays are static 
and created with a size and type dope vector. Linkage conventions 
require that actual parameters in subprogram calls and formal parameters 
in the routine called must agree in type and dimension. This precludes 
problem-dependent memory management as an automatic feature. Since array 
bounds must be predeclared, their alteration, for any reason, requires 
complete recompilation. 

Also, compilation itself is a rather slow process. The 
nominal compilation rate for APPLE II UCSD Pascal is two hundred source 
Statements per minute. This becomes tedious for the long programs devel- 
oped for this package, some of which are on the order of three thousand 
source lines (exclusive of comments). 

The final important criticism is the lack of a predefined 
exponentiation operator, a crucial primitive for nonlinear optimization 
work. For development purposes, an exponentiation operator has been 
written in UCSD Pascal. A much better solution would be a fast machine 
language function, or utilization of special arithmetic processing 


hardware. 
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E. MEMORY MANAGEMENT 

A major challenge to the microcomputer programmer attempting to code 
a complicated algorithm is memory management. The interaction between 
the two types of memory (fast "core" and slower peripherial) is a deli- 
cate affair and an eeerenely important influence on overall program 
efficiency. As a general rule, optimization programs will always expand 
to fill available core memory. Memory not occupied by program code will 
contain the data structure representation. In order to extend problem 
Size further, either the data or the program code (perhaps both) must be 
partitioned in some manner and the resulting pieces moved in and out of 
fast memory as needed. Both these approaches are feasible with today's 
microcomputers, but the intrinsic code management features of UCSD Pascal 
make program code segmenting much easier to implement. 

Data partitioning (so-called "“in-core/out-of-core" operations) 
requires additional control logic and data structures to work effectively. 
As the purpose of this research is demonstrative in nature, data parti- 
tioning experimentation has been deferred for future package enhancement. 
Only program code segmentation has been utilized to dynamically manage 
memory resources. Such a process is termed "swapping" or "overlaying." 

Due to the iterative nature of mathematical programming algorithms 
and the time expense involved in accessing current disk drives, one must 
be extremely mindful of how the program is segmented. It is desirable to 
have minimal program code in residence at all times during the critical 
solution process. If the divisions are too fine, then execution speed 
suffers as the disk is constantly being accessed. Keeping too much of 


the program in core reduces disk access at the expense of data storage. 
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Therefore, each program must be examined carefully and partitioned 


to yield the best compromise of speed and storage management. 


F. THE USER INTERFACE 
All computer programs require a certain level of human interaction 
to produce results. Depending on the particular application, this 
participation may range from minimal operations such as data input and 
output interpretation to complex control decisions. Whatever the degree 
of interaction necessary, software must be designed with the user in mind 
if maximum benefit is to be derived from its use. Programs that properly 
take the human element into consideration and attempt to promote meaning- 
ful man-machine dialog are termed "user-friendly." As microcomputers are 
interactive devices usually operated by a single person with the entire 
computing system at arm's reach, user involvement is inescapable. 
Software designed for the small computer must therefore be user-friendly. 
1. Data Input 

Thenson [Ref. 16] states that programs should be designed to 
minimize the probability of errors in the user input process. Likely 
Causes Of these errors are incompatible input format, invalid characters 
in the input field, and inadmissible input values. 

Input formats should be flexible whenever possible, ailowing 
user control to effect format reconfiguration. Otherwise, as a default, 
the easiest method for the user should prevail and internal conversion 
Can then be undertaken as needed. 

The frustrating problem of invalid characters in the input stream 
can be avoided by the use of buffering and filtering techniques to 
isolate undesired values. This error is often encountered during 
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keyboard input where, for example, an alphabetic character is entered 
when a numeric digit is required. Most operating systems respond to this 
difficulty by terminating program execution, an intolerable situation. 
Morgenson [Ref. 10] describes a procedure in detail where all numeric 
values are input as string variables. A lexicographic scan on the string 
is then performed to identify offending characters and construct the 
number for program use. We use a modification of his technique. 

Inadmissible values can be isolated by screening all input 
prior to its use. Preferably, this is performed as the data is received 
by the program. Simple precautions such as confirming that the data 
conforms to required numerical ranges (range screening) and sign restric- 
tions can prevent unexplained abnormal program terminations. 

Every attempt should be made to not only detect errors, but also 
to inform the user of their presence. This involves either providing 
facilities to correct the error and continue if possible, or executing a 
graceful program exit in the event of fatal errors. Meaningful error 
messages must be displayed in order that proper corrective action can be 
taken. 

2. Program Control 

User-friendly techniques can also be easily applied to interactive 
control during program execution. A popular method to do this involves 
menu-driven selections where various choices are displayed along with 
Simple commands. The user selects a command thereby activating either 
the desired program segment or additional sub-menus as appropriate. By 


using single keystroke commands, combined with range screening and 
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suitable error messages, avery large number of alternatives can quickly 
be considered. 

Displays should be designed so that the user has enough informa- 
tion to make decisions without being overburdened visually by extraneous 
material. Selections should be arranged with the most frequently used 
choices the easiest to invoke. A good example of this is providing 
access to default parameters only on demand. These values are then 
transparent to the user, yet readily available if changes are to be 
made. 

3. Programming Tactics 

The addition of user-friendly features invariably results in 
more complicated code and longer development time. This is a small price 
to pay in view of the utility gained. Programs that are difficult to use 
typically go unused and often must be changed, at far greater cost than 
an equivalent original design, to instill user confidence. If user- 
friendly facilities are incorporated at the conceptual phases of software 
development and sensibly blended into the program structure, their cost 
can be minimized. 

For optimization programs, there is a clear delineation between 
interface and solution modules. Even with time-shared, or dedicated 
microcomputer systems, the code supporting user-friendly facilities can 
be easily isolated from the iterative portions of the program. In this 
manner, the programmer can take advantage of these indispensable input/ 
output techniques without impairing the solution process. With proper 
memory management discipline, the only penalty is increased secondary 


(offline) storage consumption. 
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i AeaOR LT inMs 


The solution programs chosen for inclusion in this package represent 
advanced methods for solving minimum cost network flow problems. Detailed 
descriptions of the algorithms employed are scattered among various 
references, and are not available in a collected form. This section 
outlines the fundamental ideas underlying each algorithm to provide a 
Single source document. The discussions are necessarily brief and the 


reader 1s directed to the primary references for in-depth descriptions. 


A. GNET 
GNET is an extremely efficient and elegant code that solves network 
flow problems with linear costs and bounded variables. GNET uses the 
well-known primal revised simplex algorithm specialized for networks. 
1. Primal Revised Simplex Algorithm 
Consider the following linear program: 
(LP) MIN éx 


See Axe=_ 5 


© 


Kx<u 
where c is a vector of cost coefficients and A is a matrix of technological 
coefficients. Any lower bounds on the variables, x, have been eliminated 
by transformation. The matrix A may be partitioned into two sub-matrices 
B and N to yield 
A=[B: N] 
where a matrix of linearly independent columns, a basis, is represented 
by B and N consists of the remaining columns of A. 
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The vectors ¢ and x may be similarly partitioned to form 


7 ee 
& <S (C, Cy) : 


a - 
xX = (Xp Xy) ; 
which implies that 
ae. =f- 7 
CX =CXp + CX. 
Redefining variables so that every nonbasic variable is at its lower 
bound simplifies the procedure. Utilizing the notation of Bradley, Brown 
and Graves [Ref. 3], let X, be replaced by Uy > Xp (reflection) whenever 
Xp reaches its upper bound. 
Given a feasible basis, there exists a unique solution x, such that 


Bx = b, and a basic solution where gl 2 (Xp 60). 


Any solution satisfying the constraints can be written as 





e ——Ah - - 
Ax = [B N(x, Xy) = BX, + NXy = b , 
or 
ae ai 
Xp = Bb - 8B NX 
Denoting Btn as Z yields 
Xp = pth - 2x, , 
The value of the objective function ely. expressed in terms of Xy is 
Rpeeeeesl- eT -T-1,,- 
then c x = CpB b + (Cy - CpB N)X : 


Differentiating with respect to Xn provides the rate of change of the 
objective function in response to changes in Xw: 


aie 





9c X _ =I = =-Toe-l 
aX, = Cy = ChB N . 
The vector on is the solution to 
WB = és (1) 


where w is commonly referred to as the dual solution vector. 


38 








Favorable movement of the objective function (minimization) 
from a feasible solution is indicated by 
cy, > win <0 (2) 
for nonbasic variables j at zero flow. Any nonbasic variable X 
satisfying (2) will induce the following change in the solution (assuming 
all other nonbasic variables remain fixed): 


1 


ano) ab = 3x, (3) 


Ove (4) 


As this solution satisfies the explicit constraints, enforcement of the 
bound constraints will ensure that feasibility is maintained. Updating 
the variable partition completes the procedure. 
Bradley, Brown, and Graves [Ref. 3] give the following interpre- 
tation of the revised simplex procedure: 
STEP 0: Obtain a feasible solution. 
REPEAT 
STEP 1: Priceout. Select a candidate variable to enter the 
basis that satisfies (2). 
STEP 2: Ratio. Find the greatest bound such that the incoming 
variable: 
a. does not exceed its upper bound, or 
b. drives a basic variable to its lower bound, or 
c. drives a basic variable to its upper bound. 
STEP 3: Pivot. Update the solution using (3) and (4). If 
case {a) of STEP 2 applies, reflect the candidate incoming 


variable chosen in STEP 1 and leave the basis and dual 
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solution unchanged. For case (b), change the basis 
(pivot) and find a new dual solution. Case (c) requires 
that the outgoing variable first be reflected and then a 
basis and dual solution update performed. 

UNTIL STEP 1 fails to find a favorable candidate (optimality). 

2. Network Specializations 

Recall NPP? 

(NPP) MINIMIZE f(x) 


Set. Ax = b 


where A is a node-arc incidence matrix. Restricting f (x) to be a 
linear function yields a linear network programming problem (LNP): 


(LNP) MINIMIZE  ¢é!x 


cae Ax = b 


Two well-known results that characterize the bases of LNP allow 
for an extremely efficient specialization of the primal simplex procedure 
to be applied to LNP. First, all bases for LNP are formed by a set of 
columns which correspond to a spanning tree for the graph represented 
by LNP. Additionally, any such basis matrix B can always be placed in 
triangular form by simple row and column permutations. The background 
pertaining to these results can be found in any elementary network 
programming reference (e.g., Kennington and Helgason, [Ref. 26]). 

Triangulation of the basis matrix B implies that Z, which is 


i, 


defined as B ~N and thus the solution to BZ = N , can be obtained 


3Se¢ Section I. 
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by direct solution (back substitution). Also, the dual variables w 
can be found by forward substitution of (1). 

Characterization of network bases represented as spanning trees 
for the underlying graph suggests an easy method to perform the ratio 
test and pivot steps of the simplex method. Since a spanning tree by 
definition [Ref. 26] is a connected acyclic graph, any incoming nonbasic 
variable will form a cycle with the basis tree. Changes in flow as a 
result of an introduced arc will exclusively occur in unit amounts on the 
cycle in question. All off-cycle basic arcs will be unaffected. The 
direction of flow change (decrease or increase) will depend on the 
orientation of a particular basic arc with respect to the incoming arc. 
The basic arc on the cycle that reaches its bound first is removed from 
the basis (unless the incoming arc reaches its bound first). 

Bradley, Brown, and Graves [Ref. 3] employ the following 
algorithm: 

STEP 0: Given a feasible starting solution 

REPEAT 
STEP 1: Priceout. Given the dual variables w, the priceout 
becomes Se Wes Wi for nonbasic arc k that is 
directed from node i to node j. Select an incoming arc 


k with c, - W. + Wa < 0 or terminate with the 


k 
current solution as optimal. 
STEP 2: Ratio Test. Three cases analgous to bounded variable 


simplex for which: 
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a. The incoming arc reaches its opposite bound Up. 
Otherwise a basic variable is selected as the outgoing 
variable (for basic arcs on cycle with incoming arc with 
opposite orientation). 


b. Incoming arc drives basic variable down to its 


lower bound, (for basic arcs on cycle with incoming 





arc with same orientation). 
c. Incoming arc drives basic variable up to its opposite 
bound. 
| STEP 3: Update (depending upon ratio test case result). 
| a. Reflect incoming arc. 
b. Pivot update (see below). 
) c. Reflect outgoing arc, pivot update. 
Pivot Update: simultaneously perform a one-pass 
update to: 
1) modify flow on arcs in cycle by constant equal 
to ratio result, 
2) update basis tree representation (“rehang"), 
3) update dual variables for nodes whose precedessor 
chain to the root is changed by pivot operation 
(these are the rehung nodes). 
UNTIL OPTIMAL. 
The algorithm requires only integer arithmetic with most operations 
involving only addition and subtraction. Use of programming tactics and 
coding efficiencies described in Reference 3 results in an extremely 


efficient and compact code, ideal for microcomputer applications. 
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B. NLPNET 
Minimum cost nonlinear network flow problems are solved utilizing 
NLPNET [Ref. 19], a primal approach based on controlled truncation of a 
conjugate-gradient method for solving the Newton equations. As such, 
NLPNET is an extension of the unconstrained Truncated Newton methods 
described in [Ref. 27]. Operating on a maximal basis in the manner of 
Dembo and Klincewicz's Scaled Reduced Gradient method [Ref. 28], NLPNET is 
designed to take advantage of the highly desirable local convergence 
properties afforded by Newton methods while avoiding the global conver- 
gence problems and computational overhead traditionally associated with 
these procedures. These methods belong to the class of “inexact Newton 
Methods" [Ref. 29]. 
1. Basic Properties of Solutions and Methods for Nonlinear 
Optimization 
Consider the unconstrained nonlinear programming problem (NLP). 


(NLP) MINIMIZE f(x) 


xX ¢ RM 
where 

f: R"4R is a nonlinear function with the following properties. 

a. f 1S continuous and twice differentiable. 

b. At a relative minimum x*, the gradient g(x*) vanishes and the 
hessian matrix H(x*) is positive definite. 

c. For all xX eR", the level sets 

L(xo) = [x : f(x)< f(x.)] are bounded 

(e.g., Dembo and Steihaug, [Ref. 27]). 
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Algorithms to solve NLP vary considerably in specific approach, 
however, most are iterative descent methods. (Iterative means that the 
algorithms generate a series of points based on the preceeding points 
while descent implies that each new point produced by the algorithm 
improves the solution by reducing the value of the objective function. ) 

In this manner, such methods would in the ideal case converge to a 
solution point of NLP. 

For NLP, we must distinguish between two types of solution points: 
local minimum points and global minimum points. 

Definition [Ref. 30]. A point x* eR” is a local minimum point of 
NLP if there is some ce > O such that f(x) > f(x*) for all x e R" within 
distance e of x*. 

Definition [Ref. 30]. A point x* e Ruans ea global minimum point of 
NLP if f(x) > f(X*) for all X eR. 

Similarly, the concept of convergence can be viewed in a global 
and local context. Global convergence deals with the actual determina- 
tion of a solution point from an arbitrary starting point and is not to be 
confused with convergence to a global optimum. 

Definition [Ref. 30]. An algorithm is said to be globally 
convergent if for arbitrary starting points, it is guaranteed to generate 
a series of points converging to a local solution point (local minimum). 

Local convergence, on the other hand, refers to the properties of 
the algorithm in the vicinity of a solution. For example, asymptotic 


speed of convergence. One of the more important measures of this speed 


is order of convergence. 
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Definition |Ref. 30]. Let the sequence (r,} converge to r*. 
The order of convergence of (ry) is defined as the supremum of the 
non-negative numbers p satisfying 


I ear 7 TK | 


O< lim < @o, (5) 


>o - r* p 
k | ry moun 


Larger values of p in (5) imply faster convergence since the distance 


th power in a single step. 


from the limit r* is reduced by the p 
For example, quadratic convergence (p = 2) doubles precision at each 
step while for linear convergence (p = 1) the reduction ratio remains 
constant. 


Finally, if a sequence has order p convergence, then as k becomes 


very large, we have 
- rx = -~ r* p 
| co il | = 6 | Ge 0 | 
where 8 iS a constant termed the convergence ratio [Ref. 30]. 


2. Newton's Method for Unconstrained Optimization 


NLP can be solved using the well-known Newton method 
L[e.g., Ref. 30], where the function of f being minimized is approximated 
locally by a quadratic function and this approximating function is then 
minimized exactly. Thus near the point x,, we Can approximate f by the 


truncated Taylor series: 


TeX) = HeKG) + G(xo) | (x = oy) + 2 x = Xo) H(Xo) (x - Xo) 
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Minimizing the approximation by setting the derivative equal to zero 


yields: 


where g(x,) is defined as the gradient vector of the function f 
evaluated at the point x, and g'(x,) is equivalent to H(x,), the 
hessian matrix of f evaluated at x,. 

Making the appropriate symbol replacements and defining the 


vector (x - X,) aS p gives us: 
H(X,) Pp = -G(Xo) . (6) 


Equation (6) is the classical representation of the Newton method in 
an unconstrained setting with p often referred to as the Newton 
direction or the Newton step. This suggests the following iterative 
procedure given an initial guess x,: 
Step 0. k= 0. 
REPEAT 
Step 1. Solve H(x)) py = -9(x,) for Py ; 
Seep’. Xvad + Xp + Py 


UNTIL convergence. 

An equally familiar result is that under the regularity assump- 
tions given earlier this method, although well defined near a solution 
point, may not converge when far from a solution. One reason for such 
behavior is that the quadratic function is a poor approximation of f at 
the point Xr Since in its pure form the Newton direction contains 


both direction and steplength (the magnitude of the direction vector), 
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information based on the assumption of quadratic properties for f, either 
or both of these quantities may be inappropriate for the function under 
consideration. Traditionally, the method has been modified to enhance 
its convergence properties and accommodate arbitrary functions. The most 
frequent modification is the introduction of a search parameter i as 


follows: 


Xray * Xt ALPR > (7) 


where dy is chosen to minimize f in the direction Py - For points 


where the Newton quadratic approximation is "good," A, should be 


k 
approximately one. When Ay # 1, the ray search (7) is used as 

insurance that descent with respect to f is realized on the Newton 
direction. 

Finally, Newton's method can fail if H(x, ) becomes non-positive 
definite at points far from a solution. In such cases the method may not 
yield a direction of descent and can in fact seek a relative maximum 
point. The method's characteristics are summarized in Figure 4. 

3. Truncated-Newton Methods for Unconstrained Optimization 

Since the benefits of the Newton direction are greatest in the 
immediate vicinity of a solution, where the quadratic approximating 
function best describes the actual function, an accurate determination of 
the Newton Direction appears unnecessary when far from a solution. The 
use Of iterative methods to solve the Newton equations suggests a direct 


trade-off between computational effort and Newton solution accuracy. It 
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ADVANTAGES 


Locally and quadratically convergent. 


DISADVANTAGES 


The method is not globally convergent. 

It is not defined at points x, where H(x) 
is singular or essentially singular in a 
numerical sense. 

For nonconvex problems, it does not 
necessarily generate a sequence of descent 
directions. 

An n-dimensional linear system of equations 


must be solved at each iteration. 


Fig. 4. Characteristics of the Newton Method. 


is precisely this relationship that the Truncated-Newton class of methods 
[Ref. 27] seeks to exploit. At the outset of the optimization, easily 
obtainable but relatively inaccurate approximations to the Newton Direc- 
tion are tolerated with increasing accuracy demanded as Xp approaches x*. 
In order to control such a method, a measure of required accuracy 


is needed that reflects how "far" the present value of the objective 
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function is from a solution. The currently preferred measure [Ref. 27 ] 


is the relative residual 


| | ry Rie St om || where r, is defined as 


The iterative method applied to the Newton equations is therefore 
truncated when the relative residual is "small enough." 

The basic outline of a Truncated-Newton method [Ref. 27] is 
given in Figure 5. Equations (8) and (9) are conditions guaranteeing 
"Sufficient descent" [Ref. 27]. It can be shown [Ref. 31] that if these 
conditions are satisfied for dy and if {Xp} converges to an optimal 
point x* at which H(x*) is positive definite, then there is an iteration 
index k > 0 such that 4, = 1 admissible for k > k, and {Xp} converges 
Superlinearly (order > 1) to x*. These conditions therefore operate in 
conjunction with the direction finding mechanism to produce desirable 
terminal convergence properties (when the method does in fact converge). 

Using a conjugate-gradient (CG) algorithm to iteratively 
calculate the search direction Pye Dembo and Steihaug [Ref. 27] pro- 
pose the truncated conjugate-gradient algorithm (TNCG) illustrated 
in Figure 6 to serve as the minor iteration. 

There are three different ways in which the TNCG minor iteration 


can terminate: 


49 





Given an initial solution x,: 


k + 0 
REPEAT 


! Major ! 
) Iteration ) 


Compute f(X)), g(x), H(X, ) 


Test for convergence 
IF (NOT convergence) THEN 


| Minor : 
| Iteration : 


Calculate Py such that 


| 

| Seer —_ 

| H(X,) Py 7 -9(x,) 

| Using a Truncated Newton iterative method 
| 


Calculate some dy that satisfies 


f(x, + AvBL) < F(R) + ad,9(R,)'B, 5 a €(0,1/2) (8) 





g(x, F \,Dy) Pp 2 84(x,) 'D, SpE -) (9) 


7 ALP, 


esa ot 2h 


k + ktl 


| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| ENDIF 
| 


Until] convergence 


Fig. 5. Truncated-Newton Method 
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Denoting k as the major iteration counter and F(x)), g(x) ), H(x 


) 


as Fy oF and Hy respectively 


) Minor | 
| Iteration | 

| Step 1: p, + 0 
| ii oe alt, 
| d, +f, 
| ald, «rir , 
| REPEAT ; 
: Step 2: q. + Hd. 
| IF aq, < © (djd,) THEN 
| EXIT: d, if i=0 
| Pp, ot 
| K P otherwise 
| ENDIF 
| 

aie Tes 
| Step 3 a+ Pir, is d.q. 
| Bre Pa) enon 
| ey O19: 
| Sia ere oe Fe Se) 
| BAL: Pp * Pray 
| ENDIF 
| 

.— Sie 
| Step 4 B+ Peay Caay / ror. 
: Tier * Maen * 854; 

Peat aa pee 
| Tieatien * TisPasn + 85454 
| i Fel 
| 
| UNTIL EXIT 
| 
Fig. 6. Truncated Conjugate-Gradient Algorithm 
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a. The gradient vector g, points in a direction of 
negative curvature (5, 4,4, <0). In this case, the minor iteration 
returns with Pr = -G. 5 the steepest descent direction. 

b. A direction of negative curvature is encountered in the 
CG@ iteration (4H, 4, < 0) prior to satisfying the Truncated-Newton 
termination criterion. The CG procedure is terminated and the current 
estimate Pr is used. Dembo and Steihaug [Ref. 27] show that Py is a 
direction of descent. 

c. The algorithm terminates with the Truncated-Newton 
criterion. Dembo and Steihaug [Ref. 27] show that this always occurs 
in the vicinity of a strong local minimum and thus is the determining 
factor in the rate of convergence. 

Dembo and Steihaug [Ref. 27] define {n.} as the forcing 
sequence, case c above as the Truncated-Newton termination and directions 
resulting from either cases a or c as Truncated-Newton directions. 

It can be shown [Ref. 27] that the TNCG algorithm is globally 
convergent and capable of coping with regions where the hessian is 
indefinite. The following theorem indicates that the order of convergence 
can be controlled with the proper choice of the forcing sequence. 

Theorem 1. (Dembo and Steihaug [Ref. 27]) Properties of the 
Forcing Sequence. 

Assume that the Truncated-Newton iterates {X)} converge to 


x*. Then 


a2 








ae (X13 + x* superlinearly (order > 1) if (n,} + 0. 
b. {x,} + x* with order (itt) if 
- t 
me ll go Il 
for some c > 0 andO<t<tl. 


ce (x) + x* linearly (order 1) if n, is uniformly less 
than one and bounded away from zero. 


Theorem 1 indicates that by choosing 


7 - t : 
ny =e | | Sy | | (k One ) 


for some c > 0 and 0 << t < 1, the INCG algorithm will possess 
any prescribed order of convergence (1 + t) between one and two. 
The result is an adaptive algorithm that solves NLP. When far 
far from a solution, || J. || is large as is ny and little 
effort is needed to satisfy the Truncated-Newton criterion. As 
(x) } + X*, (9, } + O which implies {n,} + 0 thereby forcing {r,} + 0 and 
Py to the Newton direction. The method therefore automatically 
incorporates an increasing amount of second-order information as the 
optimization process progresses just when such information is of 
greatest use. 

4. Truncated-Newton Methods for Linearly Constrained Optimization 

Truncated-Newton methods derived for unconstrained optimization 

can be extended to solve linearly constrained nonlinear problems of the 


form 
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(LCNLP) MINIMIZE f(x) 


toe ll 


> 
~<A 
o! 


Set. 


ee | 


] 


3 


[A 


x < 
where the function f(x) is convex and separable (i.e., 
n 


foe) = fe 


j Fi (xs) and each Fa(x3) 1s a convex function). 


1 
a. The "Reduced Problem" 
In the manner of Murtagh and Saunders [Ref. 32], partition 
the columns of A as follows: 
A=[BSN], 
where the columns of B form a basis; the columns of S correspond to super- 
basic variables (nonbasic variables whose flow is allowed to vary 
between bounds); and the columns of N to nonbasic variables with flow 
fixed at either bound and not allowed to vary. Similarly, the other 
important vectors can be partitioned. 
X = [Xp Xc al (primal solution), 
g(x) = [9(Xg) 9(x<) g(xy)] (gradient vector), 
p = [Pp Pc Py (search direction vector), 





as can the function f(x) = f (Xp Xo yXw) « 


This partition allows us to re-express the constraint set 


Moe fp 
as 


TB S NI[Xp Xe XJ = b 
a B “S “N 


B Xp +5 Xe + N Xy = b 








solving for Xp (util1zing the fact that B represents a basis and thus 


-l 


= exists), 


- i} le- -l,- 
Xp 3X p =o NX : 


Now we are able to express f completely in terms of [X< Xw and denote 


= ee-5,- 6. 


this function by f. 
. ae -lz -l.- “lwo: -: 
F(X<, Xn) = f(B “b - B SX¢ - B NX wy Xe, Xv) : 


Calling f(X6, Xy) the "reduced problem" or RNLP, we note the following: 
g z f' (Xo, Xy) = [9< Dy J 
and by the chain rule, 


2 af 9 
g == 
§ aX 


—h 


aq 
A 


0 ee = 
BX Is > 9g{8 


+ 
a 
<j h 
x 
" 


B 


= af 3 
9, =z = ow + 
N 2X N gp Xn 


Therefore the modified problem becomes 


—h 
Ge 
><I| h 
Qi Q> 
ig. 
Ww 
a 

4 
(Oo 
CO 
vw, 
CO 


(RNLP) MINIMIZE f (Xe, Xw) 


: = n 
Xoo Xy € R 
S.t. 16 < Xe S$ Uc (10) 


] & Xy Su (11) 


N N° 
As the bound constraints (10) and (11) can be handled implicitly, we 
essentially are now dealing with an unconstrained problem. 
b. The Search Direction 
A generic minimization algorithm for RNLP would calculate 
successive search directions from a feasible initial point, updating x 


until convergence. Since for iteration k of this sequence, 


aig meee “tel ~ *k 


se) 


] i 





and from previous results we know that 


Ax, = [ AX. AX ¢ AX ] : 


aXy = 0 by definition, 
big = BD - BU*Sk, - Bo NK, -B b+ B sx, +B tNK, 
(#1) (k#L) (K) (K) 
: NCE oH, je Bo Saxe 
(k+1) (k) 
It follows that 
AX -B+Sax< pls 
BK) 3k) 
AX = AX = AX = I AX ; 
ue) *(k) >(k) >(k) 
AX 0 0 
Mk) 
and 
D, = AX, = I p ; 
k k , Sek) 


Given a feasible point x,, a Newton method to solve RNLP 
would involve successive solutions of a quadratic program to compute 


search directions. 


MINIMIZE 1/2p'H(x,)p + (Xe) 3 (12) 
p 
Set. P; Ss Or Xj = (jjare>) 
P; 2 0, x = u; Co) Boueieh eee 
nS 
Denoting the matrix ; as Z and altering (12) for RNLP implies 
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MINIMIZE 1/2pZ'H(%,)ZBo + 9( Xe) 'ZDe (13) 
Ps 
teres OX: = 1. (14) 


Dee mOmex es = U. . (15) 


The feasibility constraints (14) and (15) limit admissible 
search directions for variables at bounds. As Pr is chosen to be 
zero and Dc can be handled explicitly during solution of (13), only 
the components of Pp are likly to cause trouble. Dembo and Klincewicz 
[Ref. 28] point out that the major difficulty is one of wasted computational 
effort. If after Pc is determined the induced Pp violates (14) or (15), 
a new basis must be found and a new search direction [Pp Po] calculated 
until all constraints are satisfied. They further show that if one 
operates using a basis with the greatest number of variables between 
bounds (free variables), the only nonbasic variables that cause concern, 
with respect to feasibility, are those moving away from a bound. Such a 
basis is termed a maximal basis [Ref. 28]. 

The solution to (13) when H(x,) is positive definite (recall f 


aA 


and thus f are convex) is, 


t 


Z'H(X,)Zp, = -2! 


g(X_) . 
2" H(X,)Z and 2"95(x,) are often referred to as the reduced hessian and 
and reduced gradient respectively. 

We now have a system of equations in the unknown vector p 


that provides the solution to the direction finding problem for RNLP 


oy 








and hence for its equivalent, NLP. Additionally, we can recover the 
complete p vector through the relationship, 
B = ZB . 
The Truncated-Newton Congugate-Gradient method can now be applied 
directly to the quadratic approximation of RNLP. 
5. Primal Truncated-Newton Method Specialized for Networks 
The material presented to this point is totally general and 
nothing in its development relies on the fact that the problem in question 
is a network program. With this in mind, it 1s now time to examine some 
of the special structure afforded by the network representation. 
Define A of LCNLP (and RNLP) to be a node-arc incidence matrix 
of a directed network. The result is the usual formulation of the 
Capacitated minimum cost network problem with a convex separable objective 
function. We shall label these new problems NLN and RNLN. 
Any solution technique for NLN can naturally take advantage of 
the special structure of network optimization problems, dispense with 
explicit representations of the basis inverse and perform update and 
solution operations directly on a network specialization of the basis 
representation. 
First consider the form of the reduced gradient of RNLN, 
z'g = [-S'B"' 10] (4, d¢ Sy]! 
Eee 


Ml 


Ip + Yo - 
Recognizing B "5, as the estimates of the dual variables (denote 
as w) and the solution to 


t 


Bw -d, . 
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suggests that w can be determined by solving a triangular system of 
equations since network bases B (and hence their transpose) can always be 
placed in triangular form by simple permutation. Calculation of the dual 
variables in this manner would allow easy computation of the reduced 
gradient. This is a relatively standard approach to NLN. 

The solution technique just developed serves as a framework 
in which to embed an adaptive direction finding mechanism and thus 
produce a Truncated-Newton method. With the addition of control logic 
and procedures to monitor and manipulate the variable partition, we 
have all the necessary ingredients of a complete algorithm to solve 
NLN [Ref. 19]. 

Consider NLN and the associated problem RNLN. Given an initial 


feasible solution x,, 


STEP 0: Partition x, into basic, superbasic and nonbasic sets in 
such a manner that a maximal basis is established. Likewise, 
partition g and H. 
k + 0 


STEP 1: Calculate the dual variable estimates (w) by solving 


a ‘ igs 
Bw = -9, (w= -B Ip) - 


STEP 2: Compute the reduced gradient, 


at 
OX 


- 775 = -sTg-T 


= : | Se 
peas 2 9s 


O18, 





STEP 3: Test for optimality on subspace defined by active superbasic 
variables. 


ee | 24 [| < tolerance THEN optimal on current subspace 


FIND those nonbasic variables eligible to enter superbaisc 
set. Eligibility is determined by potential for feasible 
displacement from bound and subsequent reduction in objective 


function value. 


“a 


21 gilts 2 foal = 
Compute a = -N B Gp + On Nw t Gy: 


IF (af/axy > O and (xy); = 15) OR 
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(af/axy < O and (x : us) for all j € nonbasic set 


Nj 
THEN 
STOP - OPTIMAL solution has been found. 
BESE 
Add the nonbasic variable(s) not satisfying the above 
conditions to the superbasic set. 
GopOester 1. 
ENDIF 
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STEP 4: 


Sie 5: 


Sep 6: 


Ser 7: 


Ser 8: 


Compute a feasible, improving direction (Truncated-Newton 


direction) by solving 


T 


(2'HZ). = -2"§ 


using a Conjugate-Gradient (CG) algorithm with the 
termination rule 
eee 


Once 
i 2'g. 11 


kK 
where 


force, = min {forces , || 2" (on fences, < 1. t « (0,1), 


T a 


-(Z HZ(Dc) « pI g) ) 


s 
u 


CG iteration number. 


=—Je 
Tt 


Ensure all components of Dc calculated in STEP 4 
maintain feasibility. 
Desa fie enX Yoel 
Let (Ps); Ont WolG ss O and (X5) 5 ae 
(Po), = 0 if (Po) 5 2 O and (Xo) 5 = uy. 
Calculate a search direction for the basic variables by 


solving 7 Z 
BD. 2 “SPc ; 


Find a steplength A giving "sufficient descent" 
(Goldstein-Armijo conditions). 


Update the flow, 


Xe] 7% + Pp 
IF (Xp) ; at bound then pivot and replace with free arc 
from Xo (if possible) to maintain maximal basis. 
IF (x6) 5 at bound then remove from superbasic set. 
Kee tos 
GOmomolee 1. 
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Changes of bases are performed using the pivot mechanism of 
GNET [Ref. 3]. The single variable linesearch to determine A is a 
safeguarded successive cubic approximation method [Ref. 33] modified to 
incorporate the Goldstein-Armijo conditions, (8) and (9). This par- 
ticular linesearch is not crucial to the success of the method and any 
reasonable substitute could be used; although quite complex, it 1s 
reported to be very robust [Ref. 33]. The TNCG algorithm requires that 
both gradient and hessian information be available. Finally, the algo- 
rithm requires an initial feasible solution. This is not a limitation, 
but rather an advantage as the first solution can easily be provided by 
efficient linear network programs (i.e., GNET) thereby allowing the 


nonlinear code to be streamlined and overall computational effort reduced. 


C.J ENT 

ENET solves network programming problems with "elastic" ranged 
constraints and bounded variables using a modified revised primal simplex 
method. 

Elastic constraints can be violated by incurring a linear penalty 
as opposed to the rigid or inviolate constraints of the classical (network) 
programming problem. Brown and Graves [Ref. 20] point out that such 
an elastic model is a realistic and powerful portrayal of many real 
world situations. Hence it may be advantageous to relax some constraints 
(and incur a penalty) in order to satisfy others, or improve the value of 
the objective function. This class of models, therefore, offers the 


analyst complete flexibility in the formulation of NPP. 
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1. Ranged and Bounded Model 
The traditional bounded linear (network) model (LNP)? is 


(LNP) min élx 
s.t. Ax = 5 (16) 
L< 2% <0 
Addition of upper and lower ranges on the (flow conservation) 


constraints (16), yield the following ranged and bounded model (RLNP): 


(RLNP) min ex 
Bc aR <B 
1<x<u 


By introducing additional (structural) variables y and (logical) 
variables s as follows: 

y ee A, 
and 

O0<s<b-b, 
RLNP is transformed into the equivalent equality-constrained model with 
translated ranges and bounds: 

min Cytcl 


Sates S = b= Al 


s is a vector of nonnegative slack variables one for each constraint, 
that measure deviation of the current constraint value from the appropriately 


translated upper range. 


4 Section ITI, Subsection A. 
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2. Elastic Model 
Now consider the enhancement to RLNP where the constraint (flow) 
ranges are allowed to be violated. For the purposes of this paper, the 
penalty functions will be restricted to the linear case to maintain 
piecewise linearity of the objective function. Let z and z be vectors of 


the penalty coefficients (the jin 


th 


elements define the cost per unit 
violation for the i” constraint) for the upper and lower constraint 
ranges respectively. The resulting model is: 

(ERLNP) min cy+tecl+zrt+Zza 


Sateen a+ s 2 b = Al 


IR ST) 
0<s<B-b (17) 
4,7 >0 


The vectors a and r are (logical) artificial and surplus variables 
introduced to maintain equality constraints with all nonnegative variables 
as is the custom for the simplex method. 

A mathematically equivalent model for ERLNP is 


T T 


min Cy+c 1+ z(Ay - b+ Al) (18) 


with a —— 7 
Zea ND) -- anh 


z= ( -2 if Ay<b- Al 
O otherwise . 
Z, z> 0 (for convexity) 
This 1S Known as the Lagrangian form, and illuminates the fact 


that z and z are actually bounds on the variables of the dual formulation. 
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A graphical representation of the elastic cost implication for each 


constraint of ERLNP is given by Figure 7. 


Cost 






‘Region of Zero' 
Penalty 
i a 








Fig. 7. An Elastic Ranged Constraint 


In this sense, with z = z = , each constraint can be depicted as 


in Figure 8 where the discontinuities at the upper and lower ranges 


indicate that the ranges are rigid as in (RLNP) and cannot be violated 


(i.e., an infinite penalty beyond the defined range). Flow within the 


allowable range incurs zero penalty. 


Cost 


‘Region of Zero' 


Penalty 

_——— > | 

t i 
A.y¥ + Al = Ax 

b. 2 


Fig. 8. A Rigid Ranged Constraint 
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When b, = b. = b., the range collapses to a point as shown in 


—1 We 


Figure Y and the model specializes to LNP. 


Cost 


Fig. 9. Classical Equality Constraint 


Consider the relationship between the logical variables Ss, a, and 
r in ERLNP. The slack variable Ss; measures the distance a solution is 
from the upper range of constraint 1. As given by (17), S; has an 
upper bound equal to the difference between the upper and lower ranges. 
The artificial variable a; measures the distance a particular solution 
is below the lower range while the surplus variable r, measures distance 
above the upper range (both are unbounded variables in the current 
presentation). Figure 10 shows this relationship for a given constraint. 
Sig ai, and r. are mutually exclusive in any basic solution since they 
are evidently linearly dependent columns. Additionally, a. and r. must be 
equal to zero if non-basic, while S. can be non-basic at either zero or its 
upper bound b, - b.. Thus, for any given solution, the status of s., 


ass and r, may be summarized (Figure 10): 
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1. Upper range violated by r, at cost Z. Grae cress) - 


2. No violation ( ee = 0). 


it 

7m 
‘ 

oe 

=< 

I 
© 
~ ee” 


3. Lower range violated by a; at cost z. a. (s, 


Cost 








Fig. 10. Relationship Between Logical Variables for ERLNP. 


3. Elastic Primal Revised Simplex Specialized for Networks 
All solutions to ERLNP which satisfy bounds slg S x4 us) 

are feasible since the artificial and surplus variables are unbounded. 
Therefore, a solution technique for ERLNP does not need to be partitioned 
into the usual two phases: (I) where feasibility is achieved, and (ITI) 
where optimality is obtained while maintaining feasibility. Any starting 
solution satisfying bounds 1 and & will serve to begin the solution 
process. Essentially, ENET employs the primal revised simplex technique of 
GNET with modifications in logic and data structure to accommodate 


the additional arcs (variables) implied by the elastic model. Figure ll 


Shows the logical arcs employed for each node. For any solution, all 


67 





2, 


Surp lus 


Root 
- bd.) Node 








Artificial 


Fig. ll. Logical Arcs Generated for Node i 


required information for each (flow conservation) constraint and its 
penalty contribution to the objective function is generated logically 
from: 

1. which logical arc (if any) is in the basis, and 


2. which bound the slack arc S ; assumes when non-basic. 


Brown and Graves [Ref. 20] describe the following Elastic Prima! 
Simplex Network Algorithm: 


STEP O: Select a starting solution. 
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REPEAT 


STEP 1: Price Out. Given the dual variables w, the reduced 


cost of each explicit arc is: 

ey = Ww. tw, (arc k orientated 1 to j, or non- 
basic arc k at upper bound 
orientated j to 1) 


rk 


and for logical arcs, 


tw. 

Be see for a. 

+ w. 
i 0 W. for S. 
Pea We for r. 


7 j ] 
Select an incoming arc with r. < 0 or terminate with the 


current solution OPTIMAL. 

STEP 2: Ratio Test. Determine outgoing arc as in GNET. 

STEP 3: Update. Same as GNET, substituting obvious specializa- 
tions for updates involving logical arcs. 


UNTIL OPTIMAL 


The only data required by this enhanced elastic model in addition 
to classical GNET are for each node the new penalties, the difference 
between the upper and lower ranges, and whether S« is reflected. 

The resulting algorithm and data structures are well-suited 
for microcomputer implementation. The moderate increase in data storage 
1S compensated by a significant model enrichment. In particular, 

1. Total supplies and total demands no longer need be equal, nor 
must supply nodes be connected to demand nodes by paths of 
sufficient capacity to insure classical feasibility. Infeasible 


problems are intrinsically diagnosed and optimally treated. 


69 





— —— —ae 


2. All problems possess feasible primal solutions in the extended 
elastic formulations, and this provides reliable bounds for dual 
solutions. 

3. All bounded problems possess optimal solutions yielding optimal 
dual solutions. 

4. Informal relaxation (e.g., Lagrangian Relaxation) methods are 
naturally accommodated in this context. 

5. Formal decomposition methods are provided much more robust 
primal and/or dual proposals. 

6. (Mixed) integer models produce strictly feasible solutions in this 
extended Lagrangian context. 

and perhaps most important: | 

7. Solutions are reliable, inexpensive, and managerially appealing. 
Finally, as we shall see, the complete Lagrangian objective function 
yields a unifying perspective of solution properties and intrinsically 


gives profound information to the analyst or algorithms using ENET. 


D. UNET 

UNET solves elastic ranged bounded linear network problems where 
certain arcs are designated as "“l-u" or binary arcs. These arcs can 
admit flow at one of two possible values--the lower or upper bound, 


hence the label "l-u". Consider first the simplest bounded "l-u" model: 
(LULNP) = min eax 


ap. Ax = 5b 
e {1, u} (integer restriction). 


/0 








By adding the now familiar elastic range framework LULNP becomes 


T eee ra 


(ELULNP) min cyt+cC a 


1 GG 
0€s eh oh 
a, r > 0 


y= (Ye, yrs ye {l, u} (integer restriction). 


Such a formulation is termed a mixed integer problem (MIP) since 
there are both fractional (with respect to flow bounds, Xe) and integral 
(again with respect to flow bounds, x) arcs present. Note that the 
cost coefficients Cc, associated with Yi may be interpreted as fixed 
charges in the sense that for admissible solutions, the cost contribution 
is either 7 or ar for each "l-u" arc, and that each "l-u" arc is 
essentially equivalent to a binary decision variable. 

Brown and Graves [Ref. 20] employ an enumeration technique to 
solve ELUNLP which exploits the elastic model structure and produces 
excellent solutions satisfying the integer restrictions with very little 
computational effort. 

Their approach is analogous to classical branch and bound [e.g., 
Ref. 1] with the following specializations and supporting observations: 

1. Any continuous relaxation of ELULNP (with integrality violated 
by one or more arcs in Me) provides a lower bound for the 

value of an optimal solution to ELULNP. 

2. Any continuous relaxation of ELULNP can be rounded to an integer 


solution (satisfying integrality for ee) with very little 


a 


effort. Further, such rounded solutions are all admissable 

(feasible in the extended elastic sense). 

3. Restrictions of ELUNLP (with one or more arcs in ay fixed 

at a bound) are all admissible (elastically feasible) and possess 

solution values no higher than the lower bound of their respec- 

tive relaxations. 
Thus, enumeration by branch and bound may be organized: 

STEP 0: Relax (free) all integrality restrictions on he 
prepare for storage of an incumbent solution, set the 
lower bound for solution value to + @, 

REPEAT 

STEP 1: Solve ELUNLP with current restrictions. 

STEP 2: Improve lower bound on solution value, if possible. 

STEP 3: Round solution (heuristically) to satisfy integer 
restrictions. 

STEP 4: Compare rounded solution with the incumbent, replace 
incumbent if possible. 

STEP 5: IF current solution value is worse than current lower 
bound, go to STEP 7. 

STEP 6: Fix a variable. Heuristically select a free member of 
Vy and fix it at a bound. Go to STEP 1. 

STEP 7: Heuristically select STEP 8 or STEP 9. 

STEP 8: Reverse a variable. Select a previously fixed variable 
and reverse it to its opposite bound, go to STEP 1. 


If no fixed variable exists, go to STEP 9. 
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STEP 9: Backtrack. Select a reversed variable and free it. Go 
wor See. 1: 
UNTIL TERMINATION. 

Note that the heuristic decision rules involve the method of 
rounding a solution and the selection criteria for candidate variables to 
fix, reverse, and free (Figure 12). Within this framework, particular 
heuristics yield a wide variety of enumeration strategies. The particular 
Strategy chosen for microcomputer use: 

1. rounds the current restricted solution in three passes, each of 

which selects variables from a class defined in terms of 9, 

where 0 < @ < .5 and ou < y, <u or T<y < al. 

Class 1: nearly integral (0 < @< .2) 
Class 2: fractional (.2 <o0< .4) 
Class 3: ambivalent (.4 < @< .5) 
The rounding heuristic sequentially exhausts variables from each 


class and rounds using a "minimal regret function," rounding away 


from the worst penalty. 


) Fix 


y; = opposite y, = bound 
bound Free 


ae (-) 


Fig. 12. Primitive Enumeration Restrictions 
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2. variables are selected for fixing by minimal global regret, 
and 
3. variables are fixed, reversed and freed via a LIFO List operating 
exclusively on the last entry in the list: 
mix: Push fixed variable on LIFO. 
Reverse: Update status of top variable on LIFO. 
Free: Pop reversed variable from LIFO, mark it as 
freed (Figure 12). 

Also note that there are only depth and value-motivated fathoming 
rules (decision rules for reversal or backtracking); feasibility plays no 
role in the enumeration except via the value of the Lagrangian objective 
punecion. 

Finiteness of this class of enumeration methods is evident 
as long as the fix/reverse/backtrack mechanism cannot yield successive 
solutions with identical restrictions. 

What 1s not immediately apparent is the remarkable effectiveness 
of these heuristics with the elastic enumeration model. Integer solutions 
and lower bounds of excellent quality are empirically produced quite 
early in the enumeration effort, permitting routine early termination of 
the search based on an optimality tolerance or on a maximum depth (per- 
missible number of fixed variables in any restriction); tuning of the 
method is easily accomplished via these two limits and the elastic 


penalties used to express the underlying model. 
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y,; = opposite y;, = bound 
bound 
Restriction 
{4V(S.)} 
Yo = Opposite Yo = opposite 
bound bound 
Y> = bound Yo = bound 


Fathom 
by Value 
or Depth 


So = "relaxed" solution, V(S.) solution values 


Fig. 13. Generic Enumeration Tree (Fixed Order) 
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IV. PACKAGE DESCRIPTION 


The mathematical programming package described here provides an 
integrated repertoire of minimum-cost network flow models for use with a 
microcomputer. Various interface and contro] modules are also described 
which reduce the user's workload and automate the solution process. 
Extensive facilities for data file management and provisions for incor- 


porating custom problem generators complete the package. 


A. OVERVIEW 
1. Package Components 

Actually a suite of separate programs coordinated by a master 
program, the package is automated wherever possible and completely 
interactive. Designed in a highly structured manner, each program is 
modular, relatively standardized and, with minor modifications, capable 
of independent operation. Routines common to more than one program or 
subject to frequent modification reside in the system library. These 
features simplify modification or deletion of existing modules, addition 
of new programs, or even the incorporation of the entire package (or some 
subset) into a larger structure. A macro view of the package is given in 
Figure 14. 

THE APPLE II microcomputer version of the package, exclusive of 
data files, spans two disk volumes (removable floppy disks). The Pascal 
operating system, system library, master program, and al! solution 
programs reside together on a volume that is always on-line. The editor 
and associated preprocessor programs are on a second volume that can 
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Fig. 14. Package Block Diagram 





be taken off-line after the desired program is loaded into main memory. 
Data files occupy additional volumes present only for read/write opera- 
tions. This partition is dictated by the size of the package code and 
has been chosen to minimize disk manipulation requirements in a two-drive 
system (the minimum practical configuration). Package operation is thus 
able to proceed with volume exchange required only for data file access 
and editing procedures. 
2. Data Files 
a. Structure 

Although several quite distinct models are supported, a 
common format is possible through the use of Pascal's facility of variant 
records. A Pascal record 1s a compound structure of arbitrary types of 
data which, when composed of types with partly identical components, is 
termed a variant record. A portion of the record is the same for al] 
occurences while the remainder (the variant part) may differ depending on 
the value of an indicator variable (also part of the record). Wirth 
[Ref. 34] gives an excellent description of such data structures and 
their employment. This allows use of a single data structure for inter- 
program transfer of information with specific portions of the package 
extracting only those items they actually need. Each program then 
converts this input data into the required internal data representation. 
Solution technique selection can thus be accomplished without user 
intervention as the data record contains enough information to determine 
problem class and a record can be accessed by any portion of the package. 


The use of the same type record for all problems provides streamlined 
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data access procedures without sacrificing the efficiency of custom 
internal data structures for each program. 

A typical data file is a collection of randomly accessable 
records, each of which has three possible structures: (1) header, (2) 
arc, and (3) node. The Pascal interpretation of this scheme is described 
in Appendix A. Although only one header record is allowed per file, any 
number of arc and node records (subject to free disk volume space) can be 
contained in a single file without regard to order. All arcs of the 
model must be explicitly represented in the data file; however, only 
those nodes with attributes not equal to default values must be included. 
Solution programs assign appropriate default values to all nodes and then 
process the input file making note of the nodes that deviate from default 
settings. A disk volume dedicated to one file can accommodate approxi- 
mately three thousand such records (APPLE II). 

Constructed by the EDITOR during file creation, the header 
record describes the problem represented by the file. Problem name, 
problem type, number of nodes, number of arcs, and access history are 
included in this record. 

Each arc record contains the arc name, source and destination 
nodes, bounds on flow, initial flow, and cost information. Depending on 
the type of arc represented, the cost data varies. For linear cost 
functions, only the cost per unit flow and the "“I-u" status are needed. 
Nonlinear cost functions require specification of the function identifica- 
tion (assigned number in the function library) and coefficient structure. 

Node records specify name, identification number in the 


model, type, flow requirements, flow range, and penalty for range violation. 
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The identification number must be unique for a given model (enforced by 
the EDITOR on file creation and update). Flow requirements determine 
the node type as follows: zero flow--transshipment, positive flow-- 
supply, and negative flow--demand. 
Internal data representations will be presented as the 

various programs are discussed. 

b. File Naming Conventions 

In order to reduce operator workload, a menu-driven data 
file selection technique is employed. Each data file name, regardless 
of the type model it represents, contains the suffix ".net." Whenever a 
file is to be accessed, a directory for the on-line volumes is displayed, 
files with this suffix are identified, and a single keystroke selects the 
desired file. This approach ensures that only files compatible with the 
package are accessed and greatly simplifies their retrieval (from the 
user's point of view). The EDITOR program appends the required suffix to 
user-designated file names upon file creation. 
3. User Friendly Features 

Menu-driven displays and single keystroke option selection 
make the package easy to use. On input, all user entries are examined 
for errors in the sense of range, type matching (numeric vs. alpha), and 
context as applicable. 

Two general menu formats are standard throughout the package. 
The first is used for selecting courses of action. In such a situation, 
one of the displayed items must be selected before the program continues. 
The option selection menu from the EDITOR program serves as an example of 


this process and is shown in Figure 15. A common decision point in 
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EDITOR PROGRAM 


OPTIONS 
HTH RTT 


A( iter an existing file 

B(rowse through file 

C(reate a new file 

R(emove a file [permanently] 

T(ransfer a file 

<ESCAPE> EDITOR program and return to MASTER program 
OPTION DESIRED 


---[ |--- 


Fig. 15. Typical Option Selection Menu 


optimization programs concerns the setting of program parameters. The 
second menu format deals with this by displaying the default values or 
choices as appropriate. Jo change a particular value, the user enters 
the menu-designated symbo! associated with the parameter in question. If 
the parameter is an ON/OFF choice, the change is made when the symbol is 
entered; for numerical quantities, the user is prompted for the new 
value. Menu updating occurs automatically until! the user indicates that 
all values are correct. Figure 16 illustrates a typical parameter 
control menu. Choices from either format may invoke further sub-menus. 
Only selections contained in the current menu are admitted: when the 
user enters a choice symbol not depicted by the menu in use, an appro- 


priate error message is issued. 
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+ re ee er. 


ELASTICNET 


OPTIONS 
iaididimiad 
C(omplete printout --> ON 
D(etailed printout --> OFF 
F(ile output --> ORR 
H( ard copy --> OFF 
I(nternal array dump = --> Ole 2 
M(odel alteration --> OFF 
Are options satisfactory ? Y(es or N(o~ --> 


Fig. 16. Typical Parameter Control Menu 


B. CONTROL AND UTILITY PROGRAM DESCRIPTIONS 
1. Master Program 

The master program orchestrates package operation by passing 
control to appropriate programs as requested by the user (either directly 
or indirectly). This is accomplished with two levels of command internal 
to the master program as shown in Figure 17. The outer command level 
allows the choices of (1) data file manipulation, (2) problem solution, 
and (3) package exit. The first and third choices transfer control to 
the EDITOR program and the Pascal operating system respectively. A 
problem solution request activates the solution command level and package 
flow proceeds as depicted in Figure 18. 

Problem solution begins with user selection of the input file. 
Once the data file has been identified, the course of action is determined 


by problem size (number of nodes, number of arcs) and type as indicated 
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COMMAND 





| | 
| | 
| | 
| | | 
| | | 
| | Package | 
| | Exit I socution | | 
: | | CONTROL ! | 
MODULE 
| | | | | 
| | | | 
| | | | 
| | | | 
| | | 
| | 
| | 
: | | | | | 
| | EDITOR | aan | | 
| | PROGRAM | ! PROGRAMS | | 
= | | 
| | | | 
» a 


Fig. 17. Outer Command Logic 
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* indicates user interaction required (otherwise automatic) 
¢ indicates elastic & "“o-u" 
+ indicates warning message issued 


Fig. 18. Solution Process Logic 
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by the file header record. Linear problem files are those with al] 
linear arc cost functions to include problems with explicit elastic 
ranged or "l-u" constraints. Nonlinear elastic or "l-u" problems are not 
Supported and are treated as ordinary nonlinear models. 

All linear problems, size permitting, are routed to ELASTICNET 
(the combination ENET/UNET program), otherwise GNET is invoked (again 
size permitting). Nonlinear problems are first examined to ensure that 
their size is commensurate with NLPNET capabilities and that all required 
cost functions are resident in the system library. If either one of 
these tests fails, the user is given the option to linearize and use 
only GNET for an approximate solution (assuming size is within GNET 
limits). Once it has been decided that NLPNET can safely be called, 
initial flow feasibility is determined by a straightforward node flow 
conservation calculation. Feasible problems go directly to NLPNET for 
solution, however, infeasible problems are temporarily linearized and a 
feasible starting point obtained by GNET prior to being sent to NLPNET. 
In the event that GNET determines that no feasible solution exists, the 
solution process 1s terminated. Error messages are issued whenever the 
outcome of the process differs from that expected for the type file in 
question. 

A variety of elementary linearization options are offered as 
noted in Figure 19. The arc cost functions are evaluated at user desig- 
nated points and combined as specified to obtain point or interval 
approximations which then serve as cost coefficients for GNET. These 
costs are discarded after use by GNET. It should be noted that this 1s 


not a piecewise linearization and is not intended to be a complete 
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1. Point gradient eval @ LCC* of bounds (t*1 + (1-t)*u). 

2. Point gradient eval @ midpoint of bounds (t = 0.5). 

3. Point gradient eval @ LCC of lowerbound and present flow. 
4. Point gradient eval @ LCC of present flow and upperbound. 
5. Interval linearization from lowerbound to specific point. 

6. Interval linearization over entire bound interval. 

7. Assign zero cost to all arcs. 

* 


LCC = Linear Convex Combination 
OPTION DESIRED 
oe eee 


Fig. 19. Linearization Options 
solution technique. GNET can support piecewise linearizations if the 
model is explicitly described by addition of the appropriate arcs. The 
user can request that all nonlinear problems be linearized, regardless of 
feasibility status, by overriding a default parameter. 

The program calls depicted in Figure 18 are automatically gener- 
ated by the solution control module of the MASTER program. Information 
such as data file names and package control instructions are passed to 
activated programs through the Pascal operating system. Once control is 
passed to a program, this information is retrieved and the subject file 
is verified as still being on-line. In the event that the volume contain- 
ing the data file has been moved off-line, the program issues an error 
message and awaits user action. 

2. EDITOR Program 
The EDITOR program allows the user to create, alter, transfer, and 


examine data files as indicated in Figure 20. Since this portion of the 
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package requires the most user participation, every attempt has been made 
to make the program as user-friendly as possible. In addition to employ- 
ing the standard interactive techniques, the EDITOR program input format 
can be reconfigured to conform to the user's needs. This redefinition of 
the input field order can be effected at any time during the editing 
process. For example, one file can be created utilizing the default 
input order, reconfiguration performed, and a second file altered with a 
completely different input field order. Instead of reordering the raw 
data, one merely changes the order in which the program expects data 
fields to be received. Additionally, fields can be omitted entirely from 
the input format and default or user-definable constant values assigned. 
Input field order for arcs and nodes may differ. 
a. File Creation 
Files can be created from text (human readable) files, 
keyboard input, or through the use of preprocessor programs. For all 
modes of file creation, the program keeps track of the filetype by typing 
the arcs as they are received (i.e., one nonlinear arc changes a linear 
file to nonlinear). Input from text files and the keyboard is examined 
for the following inconsistencies: 
1. Prohibited node numbers (nonpositive or greater than number 
of nodes. 
2. Inconsistent upper and lower limits on quantities (lower 
limit not less than or equal to upper limit). 
3. Assignment of quantity values outside of established 


limits. 
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b. File Alteration 
Records can be added, deleted and their individual data 

fields updated. Addition of records is allowed for both node and arc 
records; however, only arcs may be deleted. This restriction 1s imposed 
because no connectivity analysis is performed by the EDITOR on the 
network representation resulting from alterations. Deletion of arcs can 
isolate portions of the network, but the solution programs accommodate 
this situation effectively. Removing nodes requires an exhaustive search 
of the data file to detect arcs that must also be deleted. This capabil- 
ity 1s not presently supported by the EDITOR. 

c. File Transfer 

Data files can be transferred to the console, printer, 
another data file, a text file, or a remote computer. The data field 
format in effect at transfer governs the order of field transmission so 
the package can be linked to data bases with diverse format requirements. 
File transfers to and from other computers are limited to text files 
only, thus the ability to convert between data files and text files has 
been provided. A stand-alone file transfer program [Ref. 35] has been 
used for such communication with great success. 
3. Preprocessor Programs 

These programs are created for specific models to assist in 
data base reduction. For example, if data on a certain network exists in 
raw form such as physical measurements on the actual entities represented 
by the nodes and arcs of the model, the user can write a custom program 
to translate such information to a form suitable for package use and then 


execute the conversion directly from the EDITOR program. A catalog of 
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existing preprocessor programs 1s maintained by the EDITOR program and 
presented to the user upon request. After selecting the new data file 
name, control is passed to the desired preprocessor and returned to the 
EDITOR upon completion of preprocessing. 

Facilities are provided by the EDITOR to maintain the preprocessor 
catalog (a data file) and manually activate uncataloged preprocessors 
(which are then automatically added to the catalog). The preprocessors 
are stand-alone programs not constrained to reside on one of the package 
volumes with no (package imposed) limits on their number or individual 


$1ze. 


C. SOLUTION PROGRAMS 

Each program is partitioned into four segments as shown in Figure 21. 
The main segment (of the particular solution program) calls subordinate 
segments into memory sequentially, thereby maximizing memory available for 
data. The input, initialization, and report segments contain all proce- 
dures that communicate with either the user or other portions of the 
package. Solution segments are therefore minimal representations of 
their respective algorithms and data structures, and are devoid of any 


nonessential information. 


| | INPUT* | 
| MAIN | INITIAL IZATION* | 
| SEGMENT | SOLUTION | 
| | REPORT* | 


| 


*User-friendly structures used 


Fig. 21. Generic Solution Program Segment Map 
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The report segment automatically updates the data file with the 
new solution and displays the results of the optimization on user- 
designated peripherals. Detailed reports of the solution process may be 
provided upon user request. Subsequent program calls are controlled by 
logic contained in the report segment. 
1. GNET 
This program is constructed so that both linear and nonlinear 
problems are processed in essentially the same manner. Two minor dif- 
ferences are accommodated by a flag that keys on the problem-type field 
of the header record. First, nonlinear files sent to GNET for solution 
have the linearized arc costs placed in the initial-flow field of the arc 
data record as GNET does not use this field for linear problem input 
data. Additionally, on obtaining a feasible solution, control is passed 
to NLPNET vice the MASTER program when GNET terminates from a nonlinear 
problem. 
a. Input Requirements 
Upon activation, GNET extracts the information listed in 
Figure 22 from the input data file. Certain compacting operations take 
place, as each arc record is processed, to reduce storage requirements. 
The arc flow range is stored as one number, the upper bound, after 
translation by the lower bound. Also, the destination node list is 
maintained with all arcs having the same destination node stored in 
contiguous locations. This allows the list to be a node-length array 
instead of an arc-length array. The arc costs and source nodes are 


read directly into arc-length arrays. 
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ARCS: 
Source node number 
Destination node number 
Lowerbound on flow 
Upperbound on flow 
Cost information 
NODES: 


Net flow 


Fig. 22. GNET Input Data Requirements 


Node information, other than net flow, is ignored as GNET 
does not support ranged constraints. In the event that total supply does 
not equal total demand, the program will terminate abnormally. GNET 
operates with integer arithmetic; therefore all input values not explic- 
itly integer are truncated upon receipt. 

b. Internal Data Structure 

GNET introduces an artificial node called the "root" to 
complete the basis tree. Three node length arrays are used to represent 
the basis tree and provide a mechanism to easily traverse the tree. The 
predecessor array defines the next node on the path from a given node to 
the root. For example, predecessor (i) contains the node number of the 


predecessor of the jun 


node. The depth array stores, for each node, 

the number of nodes on the path to the root. Finally, the traversal 
array provides for each node the next node to evaluate during forward- 
substitution of the basis. This array gives a sequence of nodes which, 
combined with the respective node predecessors, define an upper triangu- 


-th 


lation of the basic arcs. The 1° element of this array is the next 


92 








= ee OOOO eeEEEEEEEEeeeeeEEeEeeeeeeeEeEEeEeEeEeEeEeEeEeEeEeEeEeeeeEeEeEeEeEeEeEeEeEeEeEeEOEOEOEOEOEOEOEOEOEOEOEeee 


SO ——— Oo” — 


node that would be visited from node i in a classical preorder® tra- 
versal. Arrays to store flow on basic arcs and the simplex dual variables 
complete the data structure. Figure 23 illustrates the basis representa- 
tion and Figure 24 describes the various arrays. 

These arrays are used by GNET to perform the simplex opera- 
tions with elementary algebraic and logic operations. A detailed 
description of this data structure and its employment is given by 
Bradley, Brown, and Graves [Ref. 3]. 

c. Solution Segment 

An all-artificial starting technique is used to begin the 
solution process. The initial basis consists of artificial arcs between 
each node and the root node. -The flow on these arcs is set equal to the 
demand or supply requirements of the respective node. Arcs are oriented 
from the root to demand nodes and to the root from supply nodes. Assign- 
ing a relatively large cost to each artificial arc ensures that a feasible 
basis will contain no artificial arcs at positive flow. Artificial arcs 
present at the conclusion of the optimization with non-zero flow are 
explicitly identified. 

2. NLPNET 
a. Input Data Requirements 

NLPNET requires the input data illustrated in Figure 25. 

The source node, destination node, arc flow bounds, and initial flow are 
read directly into arc-length arrays. Since, in a nonlinear problem, al! 


arcs potentially can have non-zero flow at optimality, the flow array 


As defined in the computer science literature [e.g., Ref. 36]. 
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Node: 1 l 2 S 4 5 6 7 
| | 
Predecessor: i .3 5 5 3 7 -5 A | 
Traversal | 4 6 l g 3 / 5 | 
Depth | 3 2 2 3 i 2 Oa) 
| | 


---> Traversal 


Fig. 23. GNET Basis Representation 
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Problem Definition (All Arc-Length Arrays Except as Noted) 


H ( ) - A node-length array of entry points into arc-length arrays 
for each head node (e.g., all arcs directed toward head 
node i are found in locations H(i),..., H(i +1) - 1 
of arc length arrays). 

T ( ) = Tail node indices (i.e., nodes which arcs are 
directed away from). 

C ( ) - Cost per unit flow. 

CP ( ) - Upper bound (capacity) on flow. (The sign bit is used 
to indicate a reflected arc.) 


Basis Representation (All Node-length Arrays) 


P ( ) = The predecessor of each node on its backpath to the root 
node. The orientation of the basic arc connecting 
node i to its predecessor is indicated by the sign 

of P(i), a negative value indicating an arc 

(i, - P(i)), and a positive value signifying an 
arc (P(i), 7) 

IT ( ) - Preorder traversal successors. IT(i) gives the 
node number of the next node to visit in preorder 
after node i. With P( ), IT( ) defines a basis 
triangulation and can be used for substitution 
solution. 

D () - Depths of the nodes in the current basis. D(i) 
gives the length of the backpath from node i to 
the root node (used in column generation). 


Solution Representation (All Node-length Arrays) 
X ( ),- Current flow of each basic arc, and cage * minus 


CPX ( ) current flow of each basic arc (i.e., if 
Pika je< Oat tow 1s X(1))- 


Fig. 24. GNET Data Structure 
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ARCS: 
Source node number 
Destination node number 
Lowerbound on flow 
Initial flow 
Upperbound on flow 
Cost function type 
Cost function coefficients 
NODES: 


Net flow 


Fig. 25. NLPNET Input Data Requirements 


must be an arc-length structure. The bounds and the node pairs 
associated with each arc are explicitly maintained, in contrast to 
the GNET scheme where non-basic arcs are coded with a single bit to 
indicate status (at upper or lower bounds). 

Nonlinear cost function definitions reside in the SYSTEM 
LIBRARY where they can easily be modified without requiring changes to 
any of the programs comprising the package. Each function is assigned 
a unique number upon inclusion in the library, so determination of a 
function's presence is a trivial matter. This is accomplished by checking 
function identification numbers against a master list located in the 
SYSTEM LIBRARY. Also required in the function definition are analytical 
expressions for the gradient and hessian. Functions can be defined using 
up to three coefficient terms; if more than three such coefficients are 
needed, then a small coding change in NLPNET is required. 

Input node information is the same as for GNET (net flow 
only for non-default nodes). As in GNET, total supply and demand 
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conservation must be achieved or the program terminates at the input 
stage. 
b. Internal Data Structure 

NLPNET also uses a root node to complete the basis. The 
basis tree 1s depicted by the familiar predecessor, depth, and traversal 
arrays which are functionally equivalent to the corresponding GNET 
Structures. A fourth array, a reverse-traversal array, allows mobility 
Opposite to that provided by the traversal array. This array is the 
inverse of the preorder traversal and 1s used for back-substitution 
during the calculation of the basic variable search direction induced by 
a superbasic direction. The basis representation is described by 
Figure 26. Arrays to maintain the variable partition, gradient vector, 
hessian diagonal vector, search direction, and dual variable estimates 
complete the data structure. These arrays allow NLPNET to perform primal 
simplex procedures directly on the basis representation in the spirit of 
GNET. A summary of the data structure is shown in Figure 27. A complete 
description of the arrays and their use is given by Dembo [Ref. 19]. 

c. User Definable Parameters 

Real arithmetic is used extensively throughout NLPNET so 
various tolerances are necessary. Additionally, the solution process 
employs several parameters that contro] the performance, convergence, and 
accuracy characteristics of the solution process. Default settings for 
these values are built into the program, however, the user may redefine 
these settings (within established limits) if desired. Figure 28 is the 


menu that the user receives when the default settings are to be modified. 
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Node i 1 2 S 4 5 6 
| | 
Predecessor* | S 4 5 5 6 6 | 
Traversal | 6 3 1 2 4 5 | 
Inverse Traversal | 3 4 2 5 6 1 | 
Depth | 3 3 Z 2 1 Oa) | 
| | 


---> Traversal Path 

*s+*> Inverse Traversal Path 

*Negative I = predecessor (J) implies that the arc is orientated 
from node I to node J. 


Fig. 26. NLPNET Basis Representation 
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Problem Definition Arrays 


) 
) 
) 
) 
) 


ae, 


Source nodes of arcs. 

Destination nodes of arcs. 

Lower bounds on arc flow. 

Upper bounds on arc flow. 

Library identification numbers of arc objective 
functions. 

Coefficients of objective function (0-3 per arc). 
Pointers to index of coefficient array contains first 
coefficient of each objective function. 

Node net flow requirements 


Variable Partition 


LB ( ) 
ES ( ) 
LN ( ) 
ISBEST( ) 


( ) 
( ) 
LRTHRD( ) 
NDEPTH( ) 


Basic arc indices. 

Superbasic arc indices. 

Nonbasic arc indices. 

Best superbasic variable to replace a given basic 
variable. 


| Basis Representation 


Traversal array (thread). 

Predecessor array. 

Inverse of traversal array (reverse thread). 
Depth array. 


: Solution Representation 


LOW 


rman. 
a fg, gag, GO, 
eee Seren eee ene 


= U0 


Arc flow. 

Arc contribution to objective function (at FLOW). 
Gradient vector. 

Hessian diagonal vector. (Hessian is a diagonal 
matrix.) 

Search direction vector. 

Dual variable estimates. 


Fig. 27. NLPNET Data Structure 
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; | 
: | 


CURRENT PARAMETER VALUES 
GENERAL TOLERANCES 


1. Alpha = 0.0001 
2. Flow = 0.0001 
3. Gradient = 0.001 
4. Price-out = 0.0001 
LINESEARCH TOLERANCES 

5 EDS = 0.001 
Gaol = 0.02 

) eet = 0.1 

8. Epsd = 0.001 

OTHER 

9. Forcing ftn = 1.0 
10. Max Iterations = 50 


Parameters OK? : Y(es or symbol to change --> 


Fig. 28. NLPNET Default Modification Menu 


d. Reports 
The program prints iteration reports and final solution 
reports without user intervention. Either of these reports may be 
selectively disabled. In any event, the input data file's initial flow 
fields will always be updated upon normal program termination. Appendix 


B describes these reports. 


3. ELASTICNET 
This program combines both ENET and UNET into one composite 
solution module that is resident in memory for the duration of the 


optimization process. In this manner, costly disk access operations 
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resulting from UNET communication with ENET during "I-u" problem solution 
are eliminated. Of course, such a scheme requires more memory resources 
for code storage than would be the case with a partitioned format; however, 
UNET 1S avery concise code and the storage reduction is minimal. 

The selection of the appropriate solution mechanism is auto- 
matically controlled by logic internal to ELASTICNET activated by the 
presence of "l-u" arcs in the input data file. Invocation of UNET 
in response to such arcs can be suppressed from the program option menu 
and only ENET utilized to process a file with "l-u" arcs (giving only a 
relaxed solution). This option allows a single file to represent both a 
free and "l-u" model. 

a. Input Data Requirements 

ELASTICNET exercises virtually the entire input data struc- 
ture and extracts the data given in Figure 29 from the data file. For 
arcs, the source node, destination node, and unit cost are read into 
arrays identical to the GNET counterparts already described. ELASTICNET, 
unlike GNET, maintains explicit representations for each arc's upper 
and lower bounds on flow. This simplifies the program and results in a 
miniscule increase in storage requirements. Finally, the "l-u" status of 
each arc is maintained in an arc-length array with zeros indicating free 
arcs and ones marking "l-u" arcs. 

Input node information consists of five numbers for each 
node: net flow at the node, upper and lower range on flow, and penalties 
for upper and lower range violation. In the absence of input information 
for ranges and penalties, default values are assigned by the input 


module. As with the other programs, all supply and demand nodes must 
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ARCS: 
Source node number 
Destination node number 
Lowerbound on flow 
Upperbound on flow 
Cost information 
"l-u" status 

NODES: 
Net flow 
Lower range on flow 
Upper range on flow 


Lower penalty for range violation 
Upper penalty for range violation 


Fig. 29. ELASTICNET Input Data Requirements 


explicitly be present in the input file. Additionally, any transshipment 
nodes for which ranges or penalties differ from default values must be 
included. The magnitude of the default settings is user definable from 
the option menu. Flow conservation (supply = demand) is not required 
for the elastic model accommodated by ELASTICNET and therefore only a 
warning message is issued when supply does not equal demand. 
b. Internal Data Structure 

ELASTICNET maintains a basis representation (Figure 30) 
very similar to the GNET structure. Additional arrays are required to 
control the UNET enumeration process and record incumbent solutions. 


These structures are described in Figures 31 and 32. 
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Problem Definition 


Arcs: 





(All Arc-Length Arrays Except as Noted) 


A node-length array of entry points by head node into 
arc lists (sign bit indicates S. reflected). 

Tail node indices (sign bit indicates fixed arc). 
Cost per unit flow. 

Lower bounds on flow. 

Upper bounds on flow (sign bit indicates reflected arc). 


(All Node-Length Arrays) 


Lower ranges. 
Upper ranges. 
Penalties for lower range violation. 
Penalties for upper range violation. 


Basis Representation (Similar to GNET) 


Predecessor array (sign bit indicates basic arc 
orientation). 

Traversal array. 

Depth array. 

Aggregated successor array (for an aggregated node, 
the cost is stored for its basic arc predecessor; 

for a disaggregated node, the number of aggregated 
successor nodes/basic arcs is stored) [Ref. 3, p. 28]. 


Solution Representation (All Node-Length Arrays) 


X ( )y- 
BUX ( ) 


ia ) - 


Arrays with current flow of each basic arc, and 
Capacity minus current flow of each basic arc 
(j.e., P(i) = j < 0, flow is X(i)). 

Dual variables. 


Figenci. ENE} Data Structure 
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Enumeration Control 


IFIX( ) - A last-in-first-out (LIFO) structure that records the 
arcs currently fixed at a bound. IFIX(i) is the arc 
number in the current enumeration of the i 
restriction. 

LGB ( ) - Maintains (3-bit) "I-u" status of each arc. On input. 
"J-u" arcs are assigned a LGB value of one while al] 
other arcs are given a value of zero. During enumer- 
ation, nonzero values of LGB indicate bounds at which 
arcs are fixed or reversed. Fixed arcs have negative 
T ( ) while free arcs have positive T ( ) values. 


Incumbent Solution Record 


IPS ( ) = The best primal solution encountered. 
IDS ( ) - The dual solution corresponding to the best primal 
solution. 


IB ( ) - Initial right hand side of constraints. 


Fig. 32. UNET Data Structure 


c. Model Alteration 
Model alteration capability has been provided that allows 
the user to adjust arc and node parameters and resolve the modified 
problem. After viewing the new solution, the final problem attributes 
may be saved to a data file, the original file updated, or another adjust- 
ment cycle performed. 
d. Reports 
The program displays (on user designated peripherals) an 
echo print of the input file and final solution reports. Automatic 
generation of these reports is the default case, however, either report 


may be disabled. Appendix B contains examples of these reports. 
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V. COMPUTATIONAL RESULTS 


This section describes the capabilities of the package with respect 
to problem size, solution speed, and versatility. As with most applied 
mathematical programming projects, availability of suitable example 
problems hindered the testing effort. Although standard test problems 
have been described in the literature, most are either very small academic 
examples or randomly-generated problems. The difficulty with using 
randomly-generated problems is that they are unstructured and therefore 
may not adequately exploit the efficiency of a programming code designed 
to solve "real-world" problems. It has been suggested [Ref. 3] that such 
random problems may even be harder to solve than naturally occurring 
formulations. Nevertheless, it has been necessary to include some 
randoml y-generated problems to provide a complete assessment of package 
Capabilities. All randomly generated problems were produced by NETGEN 
[Ref. 37] running on an IBM 370/3033 computer and subsequently transferred 
to the microcomputer using package file transfer features (EDITOR program 
and commercial data transfer program [Ref. 35]). Additionally, a few 
small academic examples were utilized. Finally, a preprocessor has been 
written to construct models which (although fictitious) more closely resemble 
real world situations than either of the two previously mentioned 


problem classes. 


A. MAXIMUM PROBLEM SIZE 
The maximum representable formulation (in terms of numbers of arcs 
and nodes) is a function of model type, the relative proportion of arcs 
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and nodes, and the microcomputer employed. All statements concerning 
problem size apply to the APPLE II microcomputer and no attempt has been 
made to translate the results to other hardware configurations. 

Although one may specify a nominal ratio of arcs to nodes and con- 
figure the programs in that manner, 1t might be necessary to alter this 
proportion to accommodate a particular model. A limitation of this package 
is that in order to effect such an alteration, both the solution program 
pertaining to the class of models in question and the MASTER program must 
be recompiled. The coding change itself is trivial (two constants in 
each program control the partition), but the compilation process is time 
consuming. One solution to this dilemma would be to maintain multiple 
copies of the package, each with different problem size capabilities, to 
represent anticipated modeling requirements. Table 1 gives the partition 
employed during the development of the package. 

The APPLE II microcomputer has approximately 39,900 words of memory 
available for program use after the Pascal operating system has been 
loaded. The combined code of the main segment of a solution program 
(always in memory) and that of its largest subordinate segment (typically 
the solution segment) will determine the memory available for data 
Storage. In the absence of data partitioning, this also defines the 
maximum size of a representable model. The memory budget for each 
solution program is shown in Table 2. 

Data storage memory requirements can be expressed as a function of 
the number of arcs and nodes represented and will, of course, differ for 


each solution program. These relationships are as follows: 
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TABLE 1 
NOMINAL PACKAGE DATA PARTITION 


Nodes Arcs 

GNET 400 2000 

ELASTICNET 300 500 

NLPNET 50 75 
TABLE 2 


PACKAGE MEMORY USAGE 


Largest 
Main Auxiliary Available 
Segment segment for Data 
GNET 0 OF folks 27311 
ELASTICNET 6565 13898 19442 
NLPNET 6804 19941 1355 


Requirements are given in BYTES (8-bits) with 39,900 
total BYTES available for program use. 


108 








1. GNET: Smeemimeeeco oN tS. 27311, 

PeemeesolCNET: 14 | A | todo) N | < 19442, 

Sy NOP MeUe come ecmrecen| Nt <el3l55, 
where | A | signifies the number of arcs and | N | the number of nodes 
represented. A summary of this information for selected combinations of 


arcs and nodes is presented in Table 3. 


B. SOLUTION TIMES 

Solution times vary with the complexity of the model under considera- 
tion. The simplest formulation, inelastic linear models (GNET), requires 
only integer arithmetic for logical comparisions and therefore produces 
the fastest solution times for a given model size. Elastic models, 
although able to take advantage of some of the efficiencies associated 
with linear models, represent a versatile but complex formulation that 
necessitates additional work to cope with the increased information 
requirements. By specifying “l-u" arcs in a model, numerous subproblems 
(each equivalent to a single ENET run) must be solved and coordinated. 
The maximum number of such enumerations (worst case) grows exponentially 


ok ae 1, where k is the number of 


with the number of "l-u" arcs ( 
“l-u" arcs), hence solution time increases proportionally. Nonlinear models 
require the most time to achieve optimality. These models employ vast 
amounts of real arithmetic (including transcendental function computations 
which are very time consuming on the APPLE II) to perform the necessary 
functional evaluations, direction finding and linesearch solutions. 


Nonlinear solution times are also highly sensitive to the form of the 


objective functions. 
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NODES 


20 


40 
90 
60 
70 
80 
90 
100 
150 
200 
250 
300 
350 
400 
450 
900 
950 
600 
650 
700 
750 
800 
850 
900 
950 
1000 


ALLOWABLE ARC/NODE PARTITIONS 


GNET 


3368 
3346 
Boas 
280) 
S273 
3256 
3263 
Spal 
3188 
3076 
2963 
Zook 
2738 
2626 
Zols 
2401 
2288 
2176 
2063 
NS eid 
1838 
1726 
OURS 
1501 
1388 
1287 
1163 


TABLE 3 


ELASTICNET 


8) 


ARCS 


St 
No ge 
1314 
Were) 
WI 
1258 
1240 
1221 
1263 
sane 
1017 
924 
Son 
738 
645 
51518 


NUPNET 


227 
223 
ZS, 
215 
201 
207 
203 
WSS, 
isis 
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A wide variety of test problems have been examined and the solution 
results are presented in Table 4. The solution times are all quite 
reasonable (although orders of magnitude slower than those obtained with 
large computers) and reflect the expected relationships between model 
classes. An example of the accuracy for a representative nonlinear 


formulation is given in Table 5. 


SeePPMOELING FLEXIBILITY 

To demonstrate the flexibility of the package, a preprocessor program 
has been written to model a three echelon product distribution network. 
In such a model, products flow from production centers to customers via 
distribution centers with storage or handling costs incurred at the 
intermediate echelon. The object of this formulation is to meet demands 
from available supplies at minimum transportation and handling cost. 
Transportation costs are modeled by arcs between components of the 
various model echelons. Handling costs can be depicted by adding an 
additional arc at each distribution center with appropriate cost function 
and flow bounds. Figure 33 shows the basic structure of such a model. 

The preprocessor requires geographical coordinates and product flow 
requirements for each location to be modeled. All the transportation arcs 
are assigned linear cost functions proportional to the great circle 
distance between the various locations. Handling arc objective functions 
are assigned by the user to reflect the desired formulation. The user 
also specifies the minimum number of customers to be linked to each 


distribution center. The preprocessor then constructs a model as follows: 
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TABLE 5 
OPTIMAL SOLUTION COMPARISON 


Optimal Solution Flow 


DEC 20/60 
Arc APPLE II (Double Precision) 
1 52.1063 a LOo3 
2 47 .8936 47 .8937 
3 Nsiseaie 0, he yge ey 
4 33.5514 S325 o Ll) 
5 26.4450 26.4448 
6 21.4486 21.4489 
/ 5.000 5.000 
8 6.60537 6.6050 
9 6.9496 Se O02 
10 14 .6698 14.6613 
ll 10 .0000 10.0000 
IZ 8.88158 8.8898 
13 0.0000 0.0000 
14 13.3946 ie), S510, 
is 13.0504 13 .0498 
16 0.33020 ORR KY) 
17 0.0000 0.0000 
18 21.1184 21.1102 
Optimal Value of 
Objective Function 1453 .420 1453.417 


Objective function components are power, linear, and hyperbolic SIN 
functions. (Problem is NLPNET3). 
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Production Distribution 
Centers Centers Customers 


Supply —>(1) —> (;) > a G)—> @) —> demand 


handling 


Supply —>(2)—_> es og 
(3) > demand 


Fig. 33. Three Echelon Distribution Model With Handling Costs 


STEP 1. Production centers are linked with each distribution center. 

STEP 2. Handling arcs are added. 

STEP 3. Each production center is connected with the required number 
of customers selected in order of proximity. 

STEP 4. Any customer not linked to at least one distribution center 
by STEP 3 is connected to the closest distribution center. 

The following types of handling cost functions can be accommodated: 

i, kinear. 

2. “Geu". 

S.  iedwcost. 

4. Nonlinear. 

The linear formulation results, of course, in a completely linear 

model. Addition of "I-u" arcs transforms the model into one where 


distribution centers are either open at full capacity (with linear 
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handling costs) or closed. Fixed cost formulations incur an additional 
charge just for opening a distribution center (fixed cost) with a linear 
handling cost (variable cost) applied to each unit of flow. Such fixed 
cost models are generated using "l-u" arcs as shown in Figure 34. 
Finally, the handling costs may be represented by various nonlinear cost 
functions (potentially different for each handling arc). Two nonlinear 
functions that immediately come to mind as appropriate in a handling cost 
Situation are quadratic functions: 
Cost = f(x) = ax? +a +¢ 
lower bound < x < upper bound, 
and hybrid linear/quadratic functions: 
Gx lowerbound < x < changeover point, 
Cost = f(x) = 9 = —_ 
ax” + bx + (c)(changeover point) 
changeover point < x < upperbound, 
a>0O (for convexity). 
x = flow 
For demonstration purposes, a series of problems were generated with 
the following structural characteristics: 
1. five production centers, 
2. five distribution centers, 
3. twenty-five customers, and 
4. aminimum of eight customers per distribution center. 
The resulting 40-node, /73-arc model was replicated using the various 
handling arc objective functions described above. Solution times for 
these test problems are given in Table 6. 
Numerous modifications of the models are possible by exercising the 
"elastic" features of ENET and UNET. In this manner, realistic enhancements 


1S 


alive fixed cost + variable cost x (upperbound-lowerbound) 
upper bound 






{0, Upperbound} 






- Variable cost 






0, upperbound-lowerbound 


cost | 
Cost | 
| 1 


lower bound upperbound 


fixed 


flow 
Fig. 34. Modeling Fixed Costs 


such as optional surplus or shortage assignment can easily be incorporated 


into the model. 
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TABLE 6 


THREE ECHELON DISTRIBUTION MODEL SOLUTION TIMES 


(40 Nodes, 73 Arcs With 5 Handling Arcs) 


Form of Handling Arc 
Objective Function 


Linear (ENET) 


One "J-u" 
Two "l-u" 
Five "1-u" 


One Fixed Cost 
Two Fixed Cost 


Five Fixed Cost 


Quadratic 


Hybrid Linear/Quadratic 
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Solution 
Time (Seconds) 


40 


120 
NPAs 
584 


UES, 
180 
458 


65 


68 


Pivots or 
(Iterations) 


92 


256 
HO 
79 


260 
387 
960 


VI. CONCLUSIONS AND RECOMMENDATIONS 


The research effort described by this paper has been quite successful. 
Microcomputer adaptations of advanced algorithms for minimum cost network 
flow problems have been shown to be not only feasible but also practical 
for moderate-sized formulations. Additionally, an integrated repertoire 
of solution methods has been presented that illustrates the usefulness of 
such packages and serves as a prototype for their implementation on 
larger computer systems. This software fills a void in the existing 
spectrum of computer techniques for network problems by providing a 
Single package to deal with a variety of diverse modeling situations. 

Economic justification for the use of such microcomputer-based 
packages requires further investigation and could easily be the exclusive 
subject of a formal study. Certainly, the microcomputer will not replace 
larger computers, but instead will serve as a useful supplement. The 
capabilities of this software package infer that the microcomputer's 
niche lies in providing quick answers to relatively small problems. In 
this sense, a desktop computational device might be more convenient to 
use than the services of a large computer center. This is especially 
true if the software 1s user-friendly and easy to use. 

There are certain situations where the microcomputer is clearly 
the only choice. Consider remote sites where access to large computers 
is limited or nonexistent. In such a scenario, the portability of the 
smaller computer provides the analyst with a powerful mathematical 


programming capability that would otherwise not be available. Indeed, 
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microcomputers have already been used in primitive locations utilizing 
rudimentary power supplies, so such a proposal is not idle conjecture. 

The finite time horizons imposed on this project necessarily limited 
the scope of the initial study and there are many enhancements that could 
be pursued. First, the algorithms presented, although quite efficient, 
could be improved and further tuned for the microcomputer environment. 
Also, the use of partitioned data structures was not investigated as a 
means to solve larger problems with available resources. Successful use 
of such procedures could extend the usefulness of the package. Finally, 
the most obvious extension would be the inclusion of additional algorithms 
to accommodate other classes of network models. For example, a generalized 
network code recently made available to us by McBride [Ref. 38] would be 
a perfect complement to the algorithms already included. 

It is hoped that the success of the work presented here will further 
stimulate the development of additional mathematical programming software 
for microcomputer use. As smaller and more capable computers are inevi- 
table, the operations research community must be prepared to exploit 


their considerable potential. 
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APPENDIX A 


PACKAGE EXTERNAL DATA STRUCTURE (DATA FILES) 


CONST 
MX NUMCOEF = Se 
MXNAMELENGTH = 9; 
Diee 
SMALLSTRING = STRING[MXNAMELENGTH }; 
COEFARRAY = ARRAY[1..MXNUMCOEF ] OF REAL: 
RECTYPE = (HEADER,ARC NODE); 
NETWORK = (LINEAR,NONLINEAR,ELASTIC MIXEDINTEGER,GENERAL ) ; 
NODETYPE = (SUPPLY ,DEMAND , TRANSSHIPMENT) ; 
NETRECORD = RECORD 
CASE ENTITY: RECTYPE OF 
HEADER: (PROBLEMNAME >: SMALLSTRING; 
PROBLEMTYPE : NETWORK; 
NUMNODES : INTEGER; 
NUMARCS : INTEGER; 
DATECREATED : SMALLSTRING; 
DATELASTUPDATE : SMALLSTRING); 
ARC: (ARCNAME : SMALLSTRING; 
SOURCENODE : INTEGER; 
DESTINATIONNODE : INTEGER; 
LOWERBOUND : REAL; 
UPPERBOUND : REAL; 
INITIALFLOW : REAL; 
CASE ARCTYPE: NETWORK OF 
LINEAR : (UNITCOST : REAL; 
ZEROUARC : BOOLEAN); 
NONLINEAR : (TYPEFTN : INTEGER; 
NUMCOEF : 1..MXNUMCOEF); 
COEF: COEFARRAY) ; 
NODE: (NODENAME : SMALLSTRING; 
NODENUMBER : INTEGER; 
NODEKIND : NODETYPE; 
NETFLOW : REAL; 
LOWERRANGE REAL: 
UPPERRANGE : REAL; 
LOWERPENALTY : REAL; 
UPPERPENALTY lee 
END 
NETFILE : FILE OF NETRECORD: 
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APPENDIX B 
TYPICAL REPORTS 


This appendix contains three examples illustrating typical reports 


generated by the various solution programs. These examples are: 


Example 1: Linearization of a small nonlinear formulation with an 
infeasible starting solution. This serves as an example of both MASTER 


program and GNET output. 


Example 2: Typical NLPNET run of a feasible (starting solution) 


nonlinear formulation. 


Example 3: Output from a smal] linear formulation solved with ELASTICNET. 
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APPLE-NET 
Version II 


wee 


DATE: 6 SEP 81 


PROCESSING DATA:NLPNET1. NET NON LINEAR NETWORK PROGRAM. 
NUMBER OF NODES = 5 AER OF ARCS = 8 

PROBLEM IS INFEASIBLE 

LINEARIZATION PERFORMED: 

7. AsSSiqnment of zero costs for all arcs. 


THETA (POINT) = 0. 00000 
SOLUTION MODULE CHAIN OF EVENTS: 
GNET--> NLPNET--> SOLUTION... 
APPLENET - GNET MODULE 
VERSION IIIA OF 28 AUG 81 
2000 ARC, 400 NODE VERSION 
PILE: DATA: NLPNET1.NET CREATED: 1-SEP-81 UPDATED: 1-SEP-81 
PROBLEM NANE IS NLPNET-1 
NOMBER OF NODES = 5 NUMBER OF ARCS = 8 
ARC LIST... 
ARC PROM TO UNIT LOWER UPPER 
NAME NODE NODE cost BOUND BOUND 
ONE 1 2 0.00 0.00 80.00 
TWO 1 0.00 0.00 60.00 
THREE 2 3 0.00 0.00 50.00 
FOUR 2 m 0.00 0.00 80200 
FIVE 2 5 0:00 0:00 50.00 
SIX 3 mM 0.00 0.00 20.00 
SEVEN 3 5 0:00 0.00 60.00 
EIGHT 5 4 0.00 0.00 60.00 
NODE LIST ... FOR THOSE NODES EXPLICITLY IN DATA FILE 
NODE NODE NET NODE 
NAME NUMBER FLOW TYPE 
SUP1 1 100.00 SUPPLY 
UNAMED 2 0:00 TRANSSHIPMENT 
ONAMED 3 0.00 TRANSSHIPMENT 
DEM1 m -30.00 DEMAND 
DE M2 5 -70.00 DEMAND 
SUPPLIES = 100 DEMAND = 100 INITIAL COST = 0.00 
PINAL SOLUTION 
7 Pivots performed 
PROM TO FLOW P(X ) 
NODE NODE 
1 2 80 0 
1 3 20 0 
2 mn 30 0 
2 5 50 0 
5, 20 0 
OPTIMAL COST = 0.00 
Example 1 
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UPDATED 1-SEP-81 
8 


CREATED 1-SEP-81 
Number of Arcs = 


5 


50 NODE VERSION 


175 ARC 
DATE: 6 SEP 81 
DATA: NLPNET1.NET 
PROBLEM NASE IS NLPNET-1 
Number of Nodes 
ARC LIST«*.. 


FILE: 


-----OBJ ECTIVE----- 
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